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Abstract 


This report presents three methods of implementing the Dryden 
power spectral density model for atmospheric turbulence . Included are the 
equations which define the three methods and computer source code 
written in Advanced Continuous Simulation Language to implement the 
equations . Time-history plots and sample statistics of simulated 

turbulence results from executing the code in a test program are also 
presented. Power spectral densities were computed for sample sequences 
of turbulence and are plotted for comparison with the Dryden spectra . 
The three model implementations were installed in a nonlinear six-degree - 
of -freedom simulation of the High Alpha Research Vehicle airplane . 
Aircraft simulation responses to turbulence generated with the three 
implementations are presented as plots . 


Introduction 

An important part of an aircraft simulation is the capability to simulate atmospheric turbulence. Such a 
capability is useful in evaluating flight control system performance and assessing control effector activity. 
The Dryden turbulence model of power spectral density as specified in MIL_STD_1797A (ref. 1) is frequently 
used as the basis for turbulence simulation. The turbulence model is also discussed in reference 2. This 
report presents source code for three different approaches to implementing the model in a six-degree-of- 
freedom aircraft simulation and presents results obtained with these implementations. 

Symbols 


Variable: 

a 



d 

e 

f 

G di {co) 


Intermediate variable defined by equation (44); used in equation (43) 

Intermediate variable for the i-th or the j-th component of turbulence ; used in equations (52) 
and (53) as defined by equations (54); used in equations (56) as defined by equations (57) 
Intermediate variable defined by equation (44) 

Reference wing span 

Mapping constant used in the Tustin transform to relate frequencies in the digital filter 
response to frequencies in the analog response for the u-, v-, and w-components 
Mapping constant used in the Tustin transform to relate frequencies in the digital filter 
response to frequencies in the analog response for the i-th component, i-p,q,r 
Intermediate variable defined by equation (44) 

Intermediate variable defined by equation (44) 

Intermediate variable defined by equation (44) 

Intermediate variable defined by equation (44) 

Power spectral density of the i-th component of turbulence, discrete model 

Transfer function used to compute the i-th component of turbulence in the z- or the 5-domain 
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H. ( e J(oTv ) Transfer function used to compute the frequency response of H t j 

Intermediate variable defined by equation (44) 

K' kj Coefficient in difference equation in the MIL STD model for calculating the i-th or the j-th 

component of turbulence ; defined by equation (57) 

L[ Scale length of i-th component of turbulence 

M r+ , M r _ Intermediate variables defined by equation (5 1 ) 

Intermediate variable defined by equation (60) 

Sample mean of the j-th sample sequence of the i-th component of turbulence 
Mean, or expected value, of g 




~ (/■) 


m„ 


N 

R - 

R J- 

R J + 

R J-r 

R J+r 

R iM 


s 

Si(Q>) 

T 

T v 

V 
X 

x i 

Y 
Z 


Number of samples in turbulence sequence 
Intermediate variable defined by equation (74) 

Intermediate variable defined by equation (72) 

Intermediate variable defined by equation (77) 

Intermediate variable defined by equation (75) 

Intermediate variable defined by equation (73) 

Intermediate variable defined by equation (76) 

Autocorrelation function of the i-th component of turbulence 
Laplace transform variable 

Two-sided power spectral density of i-component of turbulence, continuous model 
Length of turbulence sequence ( T = NT V ) 

Sample interval in digital simulation of turbulence 
Total airspeed 

Intermediate variable defined by equation (57) 

Intermediate variable used in computation of i-th component of turbulence 

Intermediate variable defined by equation (57) 
z-transform variable 


v i ’ v i (*) Random white noise sequence used as inputs to the differential equations or difference equations 
for computation of the i-th component of turbulence 
Vj(z) z-transform of Vi(k) 

) Random sequence representing the i-th component of turbulence 
%i(z) z-transform of ^i(k) 

71 Pi ( = arc cos(-l) ) 

Single-sided power spectral density of i-th component of turbulence, continuous model 
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*i 





i 


£ 


*(*.) 


Z{m,) 

CO 


Time constant of the i-th component of turbulence 

Desired root-mean-square value of the i-th component of turbulence 

Sample standard deviation of the j-th sample sequence of the i-th component of turbulence 


Standard deviation of g 


Standard deviation of the sample standard deviation of / = ^ — £{<7; j) J 


Standard deviation of the sample mean of i 




Angular frequency 

Reciprocal of time constant for the w-component of turbulence (=l/T w ) 
Intermediate variable defined by equation (44) and equation (53) 


Subscripts: 

i 

j 

u 

v 

w 

p 

Q 

r 


Denotes turbulence component (u, v, w, p, q, or r) 

Denotes turbulence component (u, v, w, p, q, or r) 

Denotes linear velocity or turbulence component along x-axis (body axis) 
Denotes linear velocity or turbulence component along y-axis (body axis) 
Denotes linear velocity or turbulence component along z-axis (body axis) 
Denotes rotational velocity or turbulence component about x-axis (body axis) 
Denotes rotational velocity or turbulence component about y-axis (body axis) 
Denotes rotational velocity or turbulence component about z-axis (body axis) 


Superscripts: 

(j) Denotes j-th sample sequence 


Operators 

E{} 


and Notation: 

Denotes expected value of argument 


f- 1 ! • } 

N 

XH 


Denotes inverse Fourier transform of argument 
Denotes summation of argument over k from 1 to N 
Denotes absolute value of argument 


Abbreviations: 

ACSL Advanced Continuous Simulation Language 

HARV High Alpha Research Vehicle 

MIL STD Military Standard 
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PSD 


power spectral density 
root -mean -square 


rms 

ACSL Test Program (GUSTMDL) 

A test program named GUSTMDL was developed to evaluate implementations of the three models prior 
to installation in the aircraft simulation. GUSTMDL contains algorithms to simulate linear (u, v, w) and 
rotational (p, q, r) velocity components of turbulence using the three models: continuous, Tustin transform, 
and Military Standard (MIL STD) algorithms. The GUSTMDL program was written using the Advanced 
Continuous Simulation Language (ACSL) to be easily implemented into the FI 8 High Alpha Research 
Vehicle (HARV) simulation (refs. 3 and 4) used by Dynamics Control Branch (DCB) researchers. 

Equations 

The algebraic, differential, and difference equations describing the Dryden turbulence model are presented 
in this section. These equations were developed from the Dryden spectral model contained in reference 1 and 
were furnished by the Government. The equations were transformed into ACSL source code and implemented 
in the GUSTMDL program. Comments were included in the program to correlate code with the equations in 
the following section by equation number. 

The equations are separated according to each of the three turbulence models. The continuous model 
equations were implemented in the ACSL DERIVATIVE BLOCK of the GUSTMDL program. The Tustin 
model and MIL STD equations were implemented in the ACSL DISCRETE BLOCK named GUST. 

Continuous Model 

Equations (1) through (17) below define the continuous implementation of the Dryden turbulence model. 
In these equations the variables | r represent the body-axis u-, v-, w-, p-, q-, and 

r-components of turbulence, respectively. The variables V u , V v , V w , V p represent the Gaussian noise 
forcing functions for the u-, v-, w-, and p-components, respectively. 

Equations (1) and (2) are the differential equations defining the w-component of the turbulence. 

X W X W X W + H’ ^ vv Ty V w (1) 

The desired root-mean-square (rms) velocity magnitude of the w-component of the turbulence is 
represented by a w . 


L =x w + r w V3i M , = x w + r w 


( 2 ) 


where 



(3) 


The v-component can be obtained from equations (1) through (2) by substituting the subscript v for the 
subscript w. The time constant for the v-component is given in equation (4). 



The differential equation defining the u-component is 


(4) 
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( 5 ) 


A 1 , 2 

Z«=-—£u + <7 u JT^rK 

* U I V 


where 


= 


v 


Similarly, the p-component equations are 


L = ~— L+<y P I— v. 


T ^ p 

L p 


Vp T v P 


where 


1.9 


G p — r 


V 


r 

p V 

J _ V^vv^vv 

^r» — 


^ 2.6 

The differential equations defining the q-component are 


Vn n c 

Xq — TT x q “f “ tvt 

q 4b u , H 4b. 


y W 


*<3 


Vn n y. 

x n -1 Cw 

4b... q 4 b* w 


or 

z Vn ^ n j. 

^ = ~4b~^ q + 4b~^ w 

^ U W ^ U W 

A time constant for equation ( 1 2) can be defined by 


r =ik 

9 ffV 


The rms value of the q-component for equation (12) is given approximately by 


( 6 ) 


(7) 

( 8 ) 
(9) 

( 10 ) 


(ID 


( 12 ) 


(13) 


<7 


Q 


I TT 

I a w 


(14) 
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Similarly, the r-component equation is 


> VtT u 7T b 

~ ~ T7 + Ti “tv 


3b. t 


3K 


' vv “"'vv 

A time constant for equation ( 1 5) can be defined by 

■z r A 

JtV 


The rms value of the r-component for equation (15) is given approximately by 


( 15 ) 


(16) 


C7 r 



vv 


(17) 


Tustin Model 

The Tustin model can be implemented with recursive difference equations to compute components of 
turbulence. The definitions of tau in equations (3), (4), (6), (9), (13), and (16) for the continuous model were 
used in the Tustin turbulence equations. Similarly the definitions of sigma in equations (8), (14), and (17) 
apply. The Gaussian noise forcing functions are the same as used in the continuous model. 

The recursive difference equation to compute the u-component of the turbulence is 


Li k ) = - 


1 ~ CblJu 
1 + Cbl?u 


&(*-!) + 


2t„ 


1 Ty 
1 + C BL r u 


K(*) + v„(*-l)] 


(18) 


where Cg L is given by 


Cbl - 


-cot 


f *T \ 

1 V 

\ 2x uj 


(19) 


The recursive difference equation to compute the w-component is 


2 ( 

Uk)=~ A 


CO" 


'BL 


(a> w + C BL ) 




(Pw ~ £bl) 

K + c bl ) 2 


U k ~ 2) 


3 ca u 


(o) w + C BL ) 


^ ^ jv.(*) + v.(k- - 1) ■ + - C BL jv»(t - 2) 


( 20 ) 


The v-component can be obtained from equation (20) by substituting the subscript v for the subscript vv. 
The recursive difference equation to compute the p-component of the turbulence is 
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i3Sp- 1)+ T^; K w+ ^-‘>] 1211 


1 T 

C BL = — cot — 

^ r„ 2 t„ 


The difference equation for is 


$,(*) — 


1+4V^/W| ?i 


M 1 + 4b w C B L q fftV 


[U*) -«*-!)] 


Cut — COt 


tfV 7rVT v 

cot — 

\J^b w Tq 


The equations for the r-component are identical to those for the q-component except that the factor 
(4 b w /n) in equations (23) and (24) is replaced by the factor (3 b w jjt). The resulting difference equation 

for is 


(l-3b w C BL Jitv\ , 

Uk) = — | r (*-l) 

r 1 + 3 b w C BL JnV? ry 


-j ^ t [&(*)-&(*-!)] 

V[Mb w C BLr lxV)j 
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( 26 ) 


C BL =—cot 

T r 

nV 


( T* A 

l v 

y2T r j 


3b, 


-cot 


w 


nVT v 

ytb w Tq j 


The simple difference equation for the derivative of each turbulence component is 

i ( k)= &)-S( k ->) 

T \ r 


The difference equation form for the derivative of the u-component is 

C 


4(*>— ?.!“fe 4(*- i) +g .. 


2r„ 


(27) 


BL 


T v 1 + 


1 + Z uCbL 

The difference equation form for the derivative of the w-component is 


K(*)-V„(/c-l)} (28) 


u*)=~ 4 ? 




^ 1 — T C ^ 

1 T w L gL 


U*"2) + 


°wJ~ZT C BL 

1 v 


(l + t w C BL ) 

x j(l + 3 t w Cg^jv w (k) — 2 sf 3 T w C BL v w (k — 1) — ^1 — ^J 3 't W C w (k — 2)J 


(29) 


MIL STD Model 


The difference equations defining the MIL STD model were obtained from reference I. Equations (30) 
through (35) were furnished by the Government and implemented to calculate the u-, v-, w-, p-, q-, and 
r-components of turbulence. The definitions of tau and sigma for the continuous model were also used in the 
MIL STD turbulence equations. Again, the Gaussian noise forcing functions arc the same ones used in the 
continuous model. 


&(*) - fl - - l) + oJ^ v„(*) 

V T uJ 


( 2T ^ 


\ T v J 


4 T, 


^v{k) - 1-— ^ £ v (*-l) + <x v — ^v v (k) 


U*)= 


( IT ^ 
1 _ Z-V. 

V T w J 


AT, 


€w(*-1) + <TwJ—Vw(*) 


U*) = 


( T \ 

l_iv 

V T P J 


2T, 


Zp( k ~ 1 ) + <fpj — - v p( k ) 


(30) 

(31) 

(32) 

(33) 
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Uk) = 1 - i k(* - 1 ) + - W* • - 0] 


£•(*) = [ i - ^W* - 1)+ ^-[lvW-4(t - D] 


Variable Definition 

Table 1 correlates the algebraic symbols used in the equations with variable names used in coding the test 
program GUSTMDL. Variables are listed for each of the three implemented turbulence models. The equation 
number column was included as a reference to the equations listed in the previous section. Comments are 
included in the GUSTMDL code to help the reader identify the equations. 


Table 1. Correlation between Algebraic Symbol and ACSL/FORTRAN Variable 


Symbol 


Turbulence Model 


Equation 


Continuous 


MIL STD 


v p (k) 

K 

K( k ) 

v w (k) 

J P (k) 

”W 

Sr(k) 

l 

1 Jk) 


FILNP 


FILNU 


FILNV 


FILNW 


TURBP 


TURBQ 


TURBR 


TURBU 


FILNP 


FILNU 


FILNV 


FILNW 


MILPK 


MILQK 


MILRK 


21, 33 


23, 34 


25, 35 


FILU 


MILUK 


18, 30 









Table 1. Concluded 


Symbol 

Turbulence Model 

Equation 

no. 

Continuous 

Tustin 

MIL STD 

& 

TURBV 



2 

&(*) 


FELV 

MILVK 

20, 31 


TURBW 



2 

W*) 


FILW 

MILWK 

20, 32 







SIGP 

SIGP 

SIGP 

8 


SIGQ 

SIGQ 

SIGQ 

14 

<* r 

SIGR 

SIGR 

SIGR 

17 


SIGU (=TURBSIG) 

SIGU (=TURBSIG) 

SIGU (=TURBSIG) 


<*v 

SIGV (=TURBSIG) 

SIGV (=TURBSIG) 

SIGV (=TURBSIG) 


<r»- 

SIGW(=TURBSIG) 

SIGW(=TURBSIG) 

SIGW (=TURB SIG) 


T , 

TAUP 

TAUP 

TAUP 

9 

T ? 

TAUQ 

TAUQ 

TAUQ 

13 

T r 

TAUR 

TAUR 

TAUR 

16 

*u 

TAU 

TAU 

TAU 

6 

T v 

TAU 

TAU 

TAU 

4 


TAU 

TAU 

TAU 

3 






K 

BWING (=32.4) 

BWING (=32.4) 

BWING (=32.4) 


K 

TURBL (=1750.) 

TURBL (=1750.) 

TURBL (=1750.) 


K 

TURBL (=1750.) 

TURBL (=1750.) 

TURBL (=1750.) 


K 

TURBL (=1750.) 

TURBL (=1750.) 

TURBL (=1750.) 


L „ 

LP 

LP 

LP 

10 

V 

VTOT 

VTOT 

VTOT 


T 

V 

TSAMP (=0.0125) 

TSAMP (=0.0125) 

TSAMP (=0.0125) 
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Code 


The ACSL/FORTRAN source code for the GUSTMDL program (file GUSTMDL.CSL) is presented in 
this section. The program has MACRO, INITIAL, DYNAMIC, and TERMINAL BLOCKS. The 
DYNAMIC BLOCK is composed of the DERIVATIVE and DISCRETE BLOCKS. The three DISCRETE 
BLOCKS are the GUST, DISCRMS, and ASC. The continuous turbulence model is implemented in the 
DERIVATIVE BLOCK, and the Tustin and MIL STD models are implemented in the DISCRETE GUST 
BLOCK. Initialization of variables and variables that need to be calculated only one time are implemented in 
the INITIAL BLOCK. The DISCRETE BLOCK DISCRMS calls a macro to define root mean square and 
mean statistics for the variables in the three turbulence models. The DISCRETE BLOCK ASC calls 
FORTRAN subroutines ASCFM and WSTATS which are used to output variables for analysis in a format 
that is compatible with tools utilized by the DCB researchers. 


PROGRAM GUSTMDL 

t 

1 ****************************************************************** 


MACRO SECTION 


***********************************+****★**********************★+****★1 


RMS CALCULATES THE RMS RESONSE OF A RANDOM VARIABLE. 


MACRO RMS ( RMSX , X , TIV ) 

MACRO REDEFINE MSX , MSXD , MSXI 
MACRO REDEFINE MSXL, EPS 

CONSTANT MSXI = 0 . 

CONSTANT EPS = l.E-22 
MSXD = TIV* (X**2 - MSX) 

MSX = INTVC ( MSXD , MSXI ) 

MSXL = MAX ( MSX , EPS ) 

RMSX = SQRT(MSXL) 

MACRO EXIT 
MACRO END 


INVERTS TIME. USES 1/MINSTP IF TIME = ZERO. 


MACRO INVERT (TIV, T,MINSTP) 
MACRO RELABEL Ll , L2 
PROCEDURAL (TIV=T) 

IF (T. EQ. 0 . ) GOTO Ll 
TIV=1 . /T 
GOTO L2 
Ll . . CONTINUE 
TIV=1 . /MINSTP 
L2 . . CONTINUE 
END ! of procedural 
MACRO EXIT 
MACRO END 



MACRO DRMS ( RMSX , MNX , X, N ) 




MACRO TO COMPUTE SAMPLE STATISTICS OF RANDOM VARIABLE 
WTB , JCY 2-24-94 


MACRO REDEFINE SUM, SUMSQ, EPS 

CONSTANT SUM = 0 . , SUMSQ = 0 . 

CONSTANT EPS = l.E-22 
SUM = SUM + X 
SUMSQ = SUMSQ + X*X 
MNX = SUM/N 

RMSX = SQRT ( MAX ( (SUMSQ - MNX*SUM) / (MAX (N-l . , EPS) ) , 0 . ) ) 
MACRO EXIT 
MACRO END 


MACRO XCORR ( XCOR , MXM , X, Y, N, NTAU ) 

****+***+*********** + *********** + i'*i'i'ic+*i'*i r **i ( * + i' i ' i ' i ' i ' it:ki ' i ' i ' i ' i ' iri 'i' i ' i ' iri ' i ' iri ' 


! MACRO TO COMPUTE SAMPLE CROSS-CORRELATION OF RANDOM VARIABLE 
! WTB, JCY 6-2-97 

I 

MACRO REDEFINE SUMXY, SUMX, SUMY, MNX, MNY, NS, XDLY , YDLY, NDLY, I 
REAL XDLY (100) 

MACRO RELABEL LDLY 

IF (N .LT. NTAU + 1) THEN 
NDLY = N 
ELSE 

NDLY = NTAU + 1 
END IF 

DO LDLY 1=1, NDLY - 1 

XDLY (NDLY+l-I) = XDLY ( NDLY - I ) 

LDLY.. CONTINUE 
! END 

XDLY ( 1 ) = X 

CONSTANT SUMXY = 0 . , SUMX = 0 . , SUMY = 0 . , NS = 0 
IF (N .GE. NTAU + 1) THEN 
NS = NS + 1 

SUMXY = SUMXY + XDLY (NDLY) *Y 
SUMX = SUMX + XDLY (NDLY) 

SUMY = SUMY + Y 
MNX = SUMX /NS 
MNY = SUMY /NS 
MXM = MNX * MNY 
XCOR = SUMXY/NS - MXM 
END IF 

MACRO EXIT 
MACRO END 






INITIAL 




* * * * 1 
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ALGORITHM IALG 

4 


C INTERVAL 

i CINT 

\ — 1 
o 


MAXTERVAL 

i MXSTP = 

0.0125 


NSTEPS 

NSTP 

1 


INTEGER 

PAGSIZ 



CONSTANT 

PAGSIZ = 

55 


CONSTANT 

TSTP 

100.0 


PI 

ACOS ( - 1 . ) 



INTEGER NUMOUT 

$ ! DETERMINES 

CONSTANT 

NUMOUT 

= 1 


CONSTANT 

SAMPS = 0 



CONSTANT 

EPSLON = 

1 .E-22 


CONSTANT 

SDNOIS 

= 1 . , 

Sc 


TSAMP 

= .012 5, 

Sc 


TURBL 

= 1750. , 

Sc 


TURBSIG 

= 5. , 

Sc 


TURBUO 

= 0. , 

Sc 


TURBXVDO 

= 0., 

Sc 


TURBXVO 

= 0., 

Sc 


TURBXWDO 

= 0., 

Sc 


TURBXWO 

= 0. 


CONSTANT 

TURBNU 

= 0., 

Sc 


TURBNV 

= 0. , 

Sc 


TURBNW 

= 0. , 

Sc 


FILNU 

= 0. , 

Sc 


FILNV 

= 0., 

Sc 


FILNW 

= 0. 


CONSTANT 

TURBRO 

= 0. , 

Sc 


TURBQO 

= 0., 

Sc 


TURBPO 

= 0. , 

Sc 


BWING 

= 37.4 


CONSTANT 

TURBXQO 

= 0 ., 

Sc 


TURBXRO 

= 0 . 



irseed should be pos, odd, integer, 8 digits 


INTEGER IRSEED 

CONSTANT IRSEED = 28545269 

GAUSI ( IRSEED ) 

SIGU = TURBSIG 
SIGV = TURBSIG 
SIGW = TURBSIG 

SIGP = 1.9 / SQRT (TURBL * BWING) * SIGW 
APPROXIMATIONS OF EXPECTED VALUES OF STD. DEV. 


OUTPUT FOR ASC! 


! same as sigv,sigw 
! eq. 8 

(P AND R COMPONENTS) 


13 



SIGQ = SQRT (pi / ( 2 . *TURBL * BWING) ) * SIGW ! eq.14 

SIGR = SQRT ( 2 . *pi / ( 3 . *TURBL * BWING)) * SIGW ! eq.17 

CONSTANT VTOT = 400. ! VTOT = VELOCITY, FT/SEC 

use turbl= 1750. for Lu, Lv, Lw 
use tau for tauu, tauv, tauvw 


TAU = TURBL/VTOT ! eq.3,4,6 

TAUQ = 4.0 * BWING /(PI * VTOT) ! eq.13 

TAUR = 3.0 * BWING /(PI * VTOT) t eq.16 

LP = SQRT (TURBL * BWING) / 2.6 ! eq.10 

TAUP = LP/VTOT i eq.13 


TURBOMEGA = VTOT/TURBL ! same as 1 . /TAU, 1 . /TAUV, 1 . /TAUW 

WP = VTOT/LP ! same as 1 . /TAUP 

i 

! CONSTANTS FOR CALCULATIONS OF CONTINUOUS W AND V ! eq . 1 , 2 

i 

TURBK2U = TURBSIG*SQRT (2 . *TURBOMEGA/TSAMP) 

K2P = SIGP * SQRT (2. * WP/TSAMP) 

TURBK1VW = 2 . * TURBOMEGA 
TURBK2VW = TURBOMEGA* TURBOMEGA 
TURBK3VW = TAU* SQRT ( 3 . ) 

TURBK4VW = TURBSIG*SQRT (TURBOMEGA* * 3 /TSAMP) 

I 

|**********************************************************************1 
* CONSTANTS FOR CALCULATIONS of TURBULENCE VIA TUSTIN TRANSFORM 
{ 

TWO P I O VTNU = SQRT (2 . *PI /TSAMP) 

i 

i**********************************************************************| 

! U- COMPONENT 

I 

! constants used in eq.18 

i 

KFIL = TURBSIG* SQRT(1./ (PI * TURBOMEGA)) 

CFIL = TURBOMEGA/ TAN (TURBOMEGA *TSAMP/2 . ) ] eq.19 

FILKl = (TURBOMEGA - CFIL) / (TURBOMEGA + CFIL) 

FILK2 = KFIL /( 1. + (CFIL/TURBOMEGA) ) 

FILK3 = SQRT ( PI * (TURBOMEGA + CFIL) ) /SDNOIS 

I 

!**********************************************************************r 

i V-COMPONENT and W- COMPONENT 

i ★**★*★***★★**★★★★★★***★★+**★***★★****+*****★**★*★★**★★*★**★★★**★*★**★* i 

i 

! constants used in eq.20 

i 

KVW = SQRT ( 3 . *TURBOMEGA/ (2 - * PI ) ) 

WNPC = TURBOMEGA + CFIL 

FILK4 = 2 . * (TURBOMEGA*TURBOMEGA - CFIL*CFIL) / (WNPC*WNPC) 
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FILK5 = (TURBOMEGA - CFIL) * (TURBOMEGA - CFIL) / (WNPC*WNPC) 

FILK6 = CFIL + TURBOMEGA /SQRT( 3. ) 

FILK7 = 2 . *TURBOMEGA/SQRT ( 3 . ) 

FILK8 = TURBOMEGA /SQRT (3 . ) - CFIL 

FILK9 = KVW*TURBSIG / (WNPC*WNPC) 

used in eq.90 

KFILD1 = TURBSIG*SQRT ( 2 . *TAU/TSAMP) 
potential alternate (filwd2, filvd2) 

KWD1 = (1. - TAU*CFIL) / (1 . + TAU*CFIL) 

KWD2 = TURBSIG*SQRT (TAU/TSAMP) *CFIL/ ( 1 . +TAU*CFIL) **2 
CTWSR3 = SQRT (3 . ) *TAU*CFIL 

★ I*-********************************************************************! 

P-COMPONENT 
used in eq.21 

PFIL = SIGP* SQRT (2./ (TSAMP * WP) ) 

PCFIL = WP / TAN ( WP *TSAMP/2 . ) ! eq.22 

PKl = (WP - PCFIL) / (WP + PCFIL) 

PK2 = PFIL / ( 1 . + ( PCFIL/WP) ) 

Q-COMPONENT and R-COMPONENT 

*******************★***★**********************************************! 

used in eq.23 

CFILQ = (1. / TAUQ) / TAN (TSAMP / (2. * TAUQ) ) ! eq. 24 

QK1 = (1. - CFILQ * TAUQ) / (1. + CFILQ * TAUQ) 

QK2 = (CFILQ / VTOT) / (1. + CFILQ * TAUQ) 

used in eq.25 

CFILR = (1. / TAUR) / TAN (TSAMP / (2. * TAUR) ) ! eq. 26 

RK1 = (1. - CFILR * TAUR) / (1. + CFILR * TAUR) 

RK2 = (CFILR / VTOT) / (1. + CFILR * TAUR) 

**********************************************************************! 

constants for MIL STD calcs 8-28-97 

used in eqs . 30 - 35 

tovtau = tsamp/tau I (for u,w,v) 

tovtau2= 2.* tovtau 

tovtau4= 4.* tovtau 

tovtaup = tsamp/taup ! (for p) 

tovtaup2= 2 . * tovtaup 

tovtauq = tsamp/tauq i ( f or q) 
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tovtaur 

= tsamp/taur 

1 ( for 

q) 

milkl = 

(1. - tovtau) 

! (for 

U,w, V 

milk2 = 

(1. - tovtau2) 



milk3 = 

sqrt ( tovtau2 ) 



milk4 = 

sqrt ( tovtau4 ) 



milkS = 

(1. - tovtaup) 

! ( for 

p) 

milk6 - 

sqrt ( tovtaup2 ) 



milk7 = 

(1. - tovtauq) 

! (for 

q) 

milk8 = 

(1. - tovtaur) 

! (for 

r) 


piov4b=pi/ ( 4 . *bwing ) 
piov3b=pi/ (3 . *bwing) 

» 

i ***★****★★*★*★★★★★****+******★★**★★+★*★*****★★*★★**★*********+*****★** i 
i 

END ! OF INITIAL I 

i **★**★★★+**★★★★★**★★***★★**★*★★*★*******★*★*****★★******+★★+★+★★+***** i 

DYNAMIC 

i ********★**★****************★*★**■*■**** + *****★***•*★*★★*★****** + **★★★★*★ t 
i ★*★***** + ★**★★★**★★*★★★★★***★★★*★*★★ + * + *★★★*★*■*■**★■*■★***★★**★★★★★ + ★★★★★ | 

DERIVATIVE 

i********^************************************************************* i 
| 

1 CONTINUOUS TURBULENCE MODEL WTB, JCY 2-24-94 through 11-30-97 

I 

i *★★★★*★★★★***★*★★*★★*★★★*★★★★★★★**+*****★★★★*******★*★★**★★*★★+***★★*★ i 

1 U - COMPONENT 

1 ***★*★★*★****★ + ***★★★★★★***★** + ★*****★*★********★★***★*** + ** + *******★★ | 
I 

! used in eq . 5 

i 

TURBUD = - TURBOMEGA * TURBU + TURBK2U * FILNU 
TURBU = INTVC (TURBUD, TURBUO ) 

I ★*★★★★★*★★*★***★*★★*★★*★***★***★******★★**★★**********************★•*■** | 

! V - COMPONENT 

i ■****★★★*******★*•*★★★★★★**★★★★★★**★★**★*★★★****★★*****■***★★*** + ** + ***★★ t 
i 

! used in eq.2 

i 

TURBXVDD = - TURBK1VW * TURBXVD - TURBK2VW * TURBXV & 

+ TURBK4VW * FILNV 

TURBXVD = INTVC (TURBXVDD, TURBXVDO) 

TURBXVD I = TURBXVD 

TURBXV = INTVC (TURBXVDI, TURBXVO ) 

TURBV = TURBXV + TURBK3VW * TURBXVDI 

TURBVD = TURBXVDI + TURBK3VW * TURBXVDD 

! *★★***★★★★***★***********★**★★*★**★********★*********★**★*★*********** I 

I W - COMPONENT 

i ★********★**★**★*******★*★*********+**★★*★**********★**★**★*******★**★ i 

! used in eq.2 

TURBXWDD = - TURBK1VW * TURBXWD - TURBK2VW * TURBXW & 

+ TURBK4VW * FILNW 
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TURBXWD = INTVC (TURBXWDD, TURBXWDO) 

TURBXWDI = TURBXWD 

TURBXW = INTVC (TURBXWDI, TURBXWO ) 

TURBW = TURBXW + TURBK3VW * TURBXWDI 

TURBWD = TURBXWDI + TURBK3VW * TURBXWDD 

1**********************************************************************J 
! ALPHA , BETA COMPONENTS ADDED 4-28-97, 5-6-97 

t**********************************************************************! 

TURBALP = TURBW /VTOT ! ALPHA COMPONENT 

TURBBET = TURBV/VTOT 1 BETA COMPONENT 

TURBALPD = TURBWD /VTOT ! ALPHA COMPONENT 

TURBBETD = TURBVD / VTOT ! BETA COMPONENT 

t**********************************************************************! 

! P , Q, R COMPONENTS ADDED 5-15-95 JCY 

j *************************************★********************************! 

i P - COMPONENT 

!******★***************************************************************] 

! used in eq.7 

TURBPD = - WP * TURBP + K2P * FILNP 
TURBP = INTVC (TURBPD , TURBPO) 

i *********************************************************** *********** i 

! Q - COMPONENT 

t**********************************************************************! 

! used in eq.12 

TURBQD = - TURBQ / TAUQ + TURBWD / (TAUQ * VTOT) 

TURBQ = INTVC (TURBQD, TURBQO ) 

i**********************************************************************! 

! R - COMPONENT 

i ********************************************************************** ! 

j 

! used in eq.15 

i 

TURBRD = - TURBR / TAUR + TURBVD / (TAUR * VTOT) 

TURBR = INTVC (TURBRD, TURBRO ) 

j**********************************************************************! 

r 

END ! OF DERIVATIVE! 

i 

t*********************************************************************M 

j 

i **************★*****★**★**********************************************! 
j 

! DISCRETE BLOCK 

t 

t ★**★★*★**★★*★★*★★★**★***★★★★*★★★★*★★**********★*****★***★★************ l 
j 

DISCRETE GUST 

i a*********************************************************************' 
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INTERVAL TSGAUS = 0.0125 

*★***************★* + ************★****★***★★*★****★******* + *★★★*★★★**** | 
GAUSSIAN RANDOM PROCESS 

********************************************************************** | 
IF (SDNOIS.LE. EPSLON) GOTO LETA1 
FILNU = GAUSS (0 SDNOIS) 

FILNV = GAUSS (0 SDNOIS) 

FILNW = GAUSS (0 SDNOIS) 

FILNP = GAUSS (0 SDNOIS) 

GOTO LETA2 
LETA1 . . CONTINUE 
FILNU = 0. 

FILNV = 0 . 

FILNW = 0. 

FILNP = 0. 

LETA2 . . CONTINUE 

********************************************************************** i 

TURBULENCE VIA BILINEAR TRANSFORM, OR TUSTIN TRANSFORM 

********************************************************************** | 

U- COMPONENT 

********************************************************************** j 

UFILK = FILNU 
IF (T .EQ. 0.) THEN 
UFILKM1 = 0 . 

GUFILKM1 = 0. 

FILUDKM1 = 0. 

FILUD2KM1 = 0 . 

END IF 

eq. 18 

GUFILK = -FILK1 *GUFILKM1 + TWOPIOVTNU*FILK2 * (UFILK + UFILKMl ) 


★ ★★★★★★★★★★■A-********************************************************** , 

DIGITAL IMPLEMENTATION OF DERIVATIVE CALCULATIONS OF U-COMPONENT 

********************************************************************** i 


eq . 27 

FILUD = (GUFILK - GUFILKM1) /TSAMP 

eq. 2 8 

FILUD2 = - KWD1 * filud2kml & 

+ (kfildl*CFIL/ (1. + tau*CFIL) ) & 

* (ufilk - ufilkml) 5-9-97 

FILUDKM1 = FILUD 
FILUD2KM1 = FILUD2 
UFILKM1 = UFILK 
GUFILKM1 = GUFILK 
FILU = GUFILK 
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+*******★********★********★**★***************** I 


V-COMPONENT 

********************************************************************** | 

VFILK = FILNV 
IF (T .EQ. 0.) THEN 
VFILKMl = 0. 

VFILKM2 = 0. 

GVFILKM1 = 0. 

GVFILKM2 = 0. 

FILVDKM1 = 0. 

FILVD2KM1 = 0. 

FILVD2KM2 = 0. 

END IF 

FILVKM1 = GVFILKMl 

eq . 20 

GVFILK = - FILK4 *GVFILKM1 - FILK5 *GVFILKM2 & 

+ TWOPIOVTNU*FILK9 * (FILK6*VFILK + FILK7* VFILKMl & 

+ FILK8*VFILKM2) 

**********************★*★*********************************************] 
DIGITAL IMPLEMENTATION OF DERIVATIVE CALCULATIONS OF V - COMPONENT 
**********************************************************************! 


eq . 27 

FILVD = (GVFILK - GVFILKMl ) /TSAMP 
eq . 2 9 

FILVD 2 = -2 . *KWD1*FILVD2KM1 & 

- KWDl * *2 *FILVD2KM2 & 

+ kwd2 * ( ( 1 . +CTWSR3 ) *VFILK & 

- (2 . *CTWSR3 ) *VFILKM1 - ( 1 . -CTWSR3 ) *VFILKM2 ) ! 5-12-97 

GVFILKM2 = GVFILKMl 
GVFILKMl = GVFILK 
VFILKM2 = VFILKMl 
VFILKMl = VFILK 
FILVD2KM2 = FILVD2KM1 
FILVD2KM1 = FILVD 2 
FILVDKM1 = FILVD 
FILV = GVFILK 

******************★*******★★**********★*******************************[ 

W-COMPONENT 

WFILK = FILNW 
IF (T .EQ. 0.) THEN 
WFILKMl = 0. 

WFILKM2 = 0 . 

GWFILKMl = 0. 

GWFILKM2 = 0. 

FILWDKM1 = 0. 

FILWD2KM1 = 0. 
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FILWD2KM2 = 0. 

END IF 

FILWKM1 = GWFILKM1 

eq.20 

GWFILK = - FILK4 *GWFILKM1 - FILK5*GWFILKM2 & 

+ TWOPIOVTNU*FILK9* (FILK6*WFILK + FILK7 * WFI LKM1 & 

+ FILK8 *WFILKM2 ) 

DIGITAL IMPLEMENTATION OF DERIVATIVE CALCULATIONS OF W-COMPONENT 

a****************************************************************^^^ 

eq.27 

FILWD = (GWFILK - GWFILKM1) /TSAMP 


eq.2 9 


FILWD2 = -2 . *KWD1 *FILWD2KM1 & 

- KWD1 * *2 *FILWD2KM2 & 

+ kwd2 * ( ( 1 . +CTWSR3 ) *WFILK & 

- (2 . *CTWSR3) *WFILKM1 - ( 1 . -CTWSR3 ) *WFILKM2 ) 15-12-97 

GWFILKM2 = GWFILKM1 
GWFILKM1 = GWFILK 
WFILKM2 = WFILKM1 
WFILKM1 = WFILK 
FILWD2KM2 = FILWD2KM1 
FILWD2KM1 = FILWD 2 
FILWDKM1 = FILWD 
FILW = GWFILK 

**********************************************************************, 

P, Q, R COMPONENT ADDED 5-15-95 JCY 

a***************************************************************^^^ 

P-COMPONENT 

»t***********t*H***********n*»**»*********i»*t***********t******ii*»i 
PFILK = FILNP 
IF (T .EQ. 0.) THEN 
PFILKM1 = 0. 

GPFILKM1 = 0. 

END IF 

GPFILK = -PK1 *GPFILKM1 + PK2* (PFILK + PFILKM1 ) ! eq. 21 

PFILKM1 = PFILK i 5-2-97 

GPFILKM1 = GPFILK i 5-2-97 

FILP = GPFILK 

**i**t4***********t*4****t*********i**H*t****H******t*******4******t| 

Q-COMPONENT 

*************************************** *******************************, 

IF (T .EQ. 0.) THEN 
FILQKM1 = 0. 

END IF 

FILQ = -QK1*FILQKM1 


+ QK2* (FILW - FILWKM1 ) ! eq. 23 
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FILQKM1 = FILQ 


********************************************************************** \ 
R- COMPONENT 

********************************************************************** i 
IF (T .EQ. 0.) THEN 
FILRKM1 = 0. 

ENDIF 

FILR = -RK1*FILRKM1 + RK2* (FILV - FILVKM1 ) ! eq. 25 

FILRKM1 = FILR 

********************************************************************** i 

ALPHA , BETA COMPONENTS ADDED 5-7-97 

****★***************★*+★**********★***********************★***********! 

FILA = FILW/VTOT ! ALPHA COMPONENT 

FILB = FILV/VTOT ! BETA COMPONENT 

FI LAD = FILWD/VTOT 

FILBD = FILVD/VTOT 

********************************************************************** ! 

ALPHA , BETA COMPONENTS ALTERNATES ADDED 6-19-97 

*****★************★**★*★★****★* + *********★************************* + ** | 

FILAD2 = FILWD2/VTOT 
FILBD2 = FILVD2/VTOT 

**********************************************************************[ 

MIL STD IMPLEMENTATION OF CALCS OF u , v, w, p , q, r-COMPONENT 

8-28-97 

**********************************************************************! 
IF (T .EQ. 0.) THEN 
MILUKM1 = 0. 

MILVKM1 = 0. 

MILWKMl = 0. 

MILPKM1 = 0 . 

MILQKM1 = 0. 

MILRKM1 = 0. 

ENDIF 


MILUK = milkl*MILUKMl + sigu*milk3*uf ilk ! eq.30 
MILVK = milk2*MILVKMl + sigv*milk4*vf ilk ! eq.31 
MILWK = milk2*MILWKMl + sigw*milk4 *wf ilk ! eq.32 
MILPK = milk5*MILPKMl + sigp*milk6*pf ilk ! eq.33 
MILQK = milk7*MILQKMl + piov4b* (milwk-milwkml ) ! eq.34 
MILRK = milk8 *MILRKM1 + piov3b* (milvk-milvkml ) ! eq.35 

Calc derivatives 9-18-97 


MILUDK = (MILUK - MILUKM1 ) /TSAMP 
MILVDK = (MILVK - MILVKM1 ) /TSAMP 
MILWDK = (MILWK - MILWKMl ) /TSAMP 

save past values 
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MILUKM1 = MILUK 
MILVKM1 = MILVK 
MILWKM1 = MILWK 
MILPKMl = MILPK 
MILQKM1 = MILQK 
MILRKM1 = MILRK 

i 

i ********************************************************************** ^ 
i 

END ! GUST DISCRETE 

j 

i ***★********★★**★★*************************************************.*** i 

i 

DISCRETE DISCRMS 

i *★***★*★***★**★★★***★★*★**★★★************★★■*■*★****★*★**■*■*★★★* + *★★★★★★★ i 

! DISCRETE TO COMPUTE RMS OF TURBULENCE 

INTERVAL TRMS = 0.0125 
SAMPS = SAMPS + 1 . 

DRMS ( DRMSTURBU , DMNTURBU , TURBU , SAMPS) 

DRMS ( DRMSTURBV , DMNTURBV, TURBV , SAMPS) 

DRMS (DRMSTURBW, DMNTURBW , TURBW , SAMPS) 

i 

DRMS (DRMS FI LNU, DMNFILNU, FILNU , SAMPS) 

DRMS (DRMSFILNV, DMNFILNV, FILNV, SAMPS) 

DRMS (DRMSFILNW, DMNFILNW, FILNW, SAMPS) 

i 

DRMS (DRMSFILNP, DMNFILNP, FILNP , SAMPS) 

i 

DRMS (DRMSFILU, DMNFILU, FILU, SAMPS) 

DRMS (DRMSFILV, DMNFILV, FILV, SAMPS) 

DRMS (DRMSFILW, DMNFILW, FILW, SAMPS) 

; 

DRMS ( DRMSTURBUD , DMNTURBUD , TURBUD , SAMPS) 

DRMS ( DRMSTURBVD , DMNTURBVD , TURBVD , SAMPS ) 

DRMS ( DRMSTURBWD , DMNTURBWD , TURBWD , SAMPS ) 

DRMS (DRMSTURBP, DMNTURBP , TURBP, SAMPS) 

DRMS ( DRMSTURBQ , DMNTURBQ , TURBQ , SAMPS ) 

DRMS (DRMSTURBR, DMNTURBR, TURBR, SAMPS) 

i 

DRMS (DRMSFILUD, DMNFILUD, FILUD , SAMPS) 

DRMS (DRMSFILVD, DMNFILVD, FILVD, SAMPS) 

DRMS (DRMSFILWD, DMNFILWD, FILWD, SAMPS) 

i 

DRMS (DRMSFILP, DMNFILP, FILP, SAMPS) 

DRMS (DRMS FI LQ, DMNFILQ, FILQ, SAMPS) 

DRMS (DRMS FI LR, DMNFILR, FILR, SAMPS) 

i 

! FOLLOWING ADDED 5-8-97 
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DRMS (DRMSTURALP , DMNTURBALP, TURBALP, SAMPS) 

DRMS ( DRMSTURBBET , DMNTURBBET , TURBBET, SAMPS) 

I 

DRMS (DRMSTURALPD , DMNTURBALPD, TURBALPD, SAMPS) 

DRMS (DRMSTURBBETD, DMNTURBBETD, TURBBETD , SAMPS) 

I 

DRMS (DRMSFILA , DMNFILA , FILA , SAMPS) 

DRMS (DRMSFILB , DMNFILB , FILB r SAMPS) 

DRMS (DRMSFI LAD, DMNFILAD, FILAD, SAMPS) 

DRMS (DRMSFILBD, DMNFILBD, FILBD, SAMPS) 

! FILAD2 , FILBD 2 ADDED 6-29-97 

DRMS (DRMSFILAD2 , DMNFILAD2 , FILAD2 , SAMPS) 

DRMS (DRMSFILBD2 , DMNFILBD2 , FILBD2, SAMPS) 

I 

DRMS ( DRMSTURBQ2 , DMNTURBQ2 , TURBQ2 , SAMPS) 

DRMS ( DRMSTURBR2 , DMNTURBR2 , TURBR2 , SAMPS) 

i 

; 

DRMS (DRMSFILUD2 , DMNFILUD2 , FILUD2 , SAMPS) 

DRMS ( DRMSFILVD2 , DMNFILVD2 , FILVD2 , SAMPS) 

DRMS (DRMSFILWD2, DMNFILWD2 , FILWD2 , SAMPS) 

I 

XCORR ( XCORTW , MXMTW , TURBW, TURBW , SAMPS, 1) 

XCORR ( XCORFW , MXMFW , FILW, FILW, SAMPS, 1) 

I 

XCORR (XCORTV, MXMTV , TURBV, TURBV, SAMPS, 1) 

XCORR ( XCORFV , MXMFV, FILV, FILV, SAMPS, 1) 

I 

XCORR (XCORTQW, MXMTQW , TURBQ , TURBW, SAMPS, 0) 

XCORR (XCORFQW, MXMFQW, FILQ, FILW, SAMPS, 0) 

I 

XCORR (XCORTRV, MXMTRV , TURBR, TURBV, SAMPS, 0) 

XCORR (XCORFRV, MXMFRV , FILR, FILV, SAMPS, 0) 

! MIL STD calcs, u , v, w, p, q, r ADDED 8-28-97 

DRMS ( DRMSMILUK , DMNMILUK, MILUK, SAMPS) 

DRMS (DRMSMILVK, DMNMILVK, MILVK, SAMPS) 

DRMS (DRMSMILWK, DMNMILWK, MILWK, SAMPS) 

DRMS (DRMSMILPK, DMNMILPK, MILPK, SAMPS) 

DRMS (DRMSMILQK, DMNMILQK, MILQK, SAMPS) 

DRMS (DRMSMILRK, DMNMILRK, MILRK, SAMPS) 

1 

DRMS (DRMSMILUDK, DMNMILUDK, MILUDK, SAMPS) 

DRMS (DRMSMILVDK, DMNMILVDK, MILVDK, SAMPS) 

DRMS (DRMSMILWDK, DMNMILWDK, MILWDK , SAMPS) 

I 

END ! of DISCRMS DISCRETE 

i ********************************************************************** 


DISCRETE TO SAVE DATA TO AN ARRAY IN DRYDEN GETDATA FORMAT 
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1 


DISCRETE ASC 

INTERVAL TSASC = 0.0125 


CALL ASCFM ( T , 


TURBU , 

TURBV, 

FILV, 

FILW, 

TURBWD , 

FI LUD , 

FILNU, 

FILNV, 

TURBP , 

TURBQ, 

FILP, 

FILQ, 

TURBBET, 

TURBALPD, 

FILB , 

FILAD, 

TURBR2 , 

FILUD2, 

FI LAD 2 , 

FILBD2 , 

MILWK, 

MILPK, 

MILUDK, 

MILVDK, 

NUMOUT , 

TSTP) 


TERMT ( T . GE . TSTP ) 

END ! ASC DISCRETE 

END ! OF DYNAMIC! 

i 

TERMINAL 


St 


TURBW , 

FILU , 

Sc 

TURBUD , 

TURBVD , 

Sc 

FILVD, 

F I LWD , 

Sc 

FILNW, 

FILNP, 

Sc 

TURBR, 


Sc 

FILR , 

TURBALP, 

Sc 

TURBBETD , 

FILA, 

St 

FILBD , 

TURBQ2 # 

St 

FILVD2 , 

FI LWD 2 , 

St 

MILUK, 

MILVK , 

Sc 

MILQK, 

MILRK, 

Sc 

MILWDK, 


Sc 


CALL WSTATS (T, & 

DRMSTURBU , DRMSTURBV , DRMSTURBW , & 

DRMSTURBP , DRMSTURBQ , DRMSTURBR , & 

DRMSTURBQ2 , DRMSTURBR2 , & 

DRMSFILU , DRMSFILV, DRMSFILW, & 

DRMSFILP , DRMSFILQ, DRMSFILR, & 

DRMSFILNU, DRMSFILNV, DRMSFILNW, DRMSFILNP, & 
DRMSTURBUD , DRMSTURBVD , DRMSTURBWD, & 

DRMSFILUD, DRMSFILVD, DRMSFILWD, & 

DRMSFILUD2 , DRMSFILVD2 , DRMSFILWD2 , & 

DRMSTURAL P , DRMSTURBBET , & 

DRMSTURALPD , DRMSTURBBETD , & 

DRMSFILA, DRMSFILB, DRMSFILAD, DRMSFILBD, & 

DRMSFILAD2 , DRMSFILBD2 , & 

DMNTURBU , DMNTURBV , DMNTURBW , & 

DMNTURB P , DMNTURBQ , DMNTURBR , & 

DMNTURBQ2 , DMNTURBR2 , & 

DMNFILU, DMNFILV, DMNFILW, & 

DMNFILP, DMNFILQ, DMNFILR, & 

DMNFILNU, DMNFILNV, DMNFILNW, DMNFILNP, & 

DMNTURBUD , DMNTURBVD , DMNTURBWD , & 

DMNF I LUD , DMNF I LVD , DMNF I LWD , & 

DMNFILUD2 , DMNF I LVD 2 , DMNF I LWD 2 , & 

DMNTURBALP , DMNTURBBET , & 

DMNTURBAL PD , DMNTURBBETD , & 

DMNF I LA, DMNF I LB, DMNF I LAD, DMNF I LBD, & 
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DMNFILAD2 , DMNFILBD2 , & 
DRMSMILUK , DMNMILUK, DRMSMILVK, DMNMILVK, & 
DRMSMILWK, DMNMILWK, DRMSMILPK, DMNMILPK, & 
DRMSMILQK, DMNMILQK, DRMSMILRK, DMNMILRK , & 


DRMSMILUDK, DMNMILUDK, DRMSMILVDK, DMNMILVDK, & 
DRMSMILWDK, DMNMILWDK) 

I 

END ! END OF TERMINAL 

i 

i 

END ! OF PROGRAM ! 

SUBROUTINE WSTATS (T, 
c 

c ******* 35 + 2 = 37 rms ******************************** 


c 


* 

DRMSTURBU , 

DRMSTURBV, 

DRMSTURBW , 


★ 

DRMSTURBP , 

DRMSTURBQ, 

DRMSTURBR , 


* 

DRMSTURBQ2 , 

DRMSTURBR2 , 



* 

DRMSFILU, 

DRMSFILV, 

DRMSFILW, 


* 

DRMSFILP, 

DRMSFILQ, 

DRMSFILR, 


* 

DRMSFILNU , 

DRMSFILNV, 

DRMSFILNW, 

DRMSFILNP, 

★ 

DRMSTURBUD, 

DRMSTURBVD, 

DRMSTURBWD , 


★ 

DRMSFILUD, 

DRMSFILVD , 

DRMSFILWD , 


★ 

DRMSFILUD2 , 

DRMSFILVD2 , 

DRMSFILWD2 , 


* 

DRMSTURALP , 

DRMSTURBBET , 



★ 

DRMSTURALPD , 

DRMSTURBBETD , 



* 

DRMSFILA, 

DRMSFILB , 

DRMSFILAD, 

DRMSFILBD, 

* 

DRMSFILAD2 , 

DRMSFILBD2, 



c 

* 35 + 2 = 37 

means + 12 mils= 49 values 


******************************** 



c 

* 

DMNTURBU , 

DMNTURBV , 

DMNTURBW, 


* 

DMNTURBP, 

DMNTURBQ , 

DMNTURBR, 


* 

DMNTURBQ2 , 

DMNTURBR2 , 



* 

DMNFILU, 

DMNFILV, 

DMNFILW, 


* 

DMNFILP, 

DMNFILQ, 

DMNFILR, 


* 

DMNFILNU , 

DMNFILNV, 

DMNFILNW, 

DMNFILNP, 

★ 

DMNTURBUD, 

DMNTURBVD , 

DMNTURBWD , 


* 

DMNFILUD, 

DMNFILVD, 

DMNFILWD, 


* 

DMNFILUD2 # 

DMNFILVD2 , 

DMNFILWD2 , 


* 

DMNTURBALP , 

DMNTURBBET , 



* 

DMNTURBALPD, 

DMNTURBBETD , 



* 

DMNFILA, 

DMNFILB , 

DMNFILAD, 

DMNFILBD, 

* 

DMNFILAD2 , 

DMNFILBD2 , 



* 

DRMSMILUK, 

DMNMILUK, 

DRMSMILVK, 

DMNMILVK, 

★ 

DRMSMILWK, 

DMNMILWK, 

DRMSMILPK, 

DMNMILPK, 

* 

DRMSMILQK, 

DMNMILQK, 

DRMSMILRK, 

DMNMILRK, 

* 

DRMSMILUDK, 

DMNMILUDK, 

DRMSMILVDK, 

DMNMILVDK, 

* 

DRMSMILWDK, 

DMNMILWDK) 




C****************************************************************** 

C* SUBROUTINE WSTAT 
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C* WRITES AN OUTPUT FILE IN THE GETDATA FORMAT ASC . to IMPORT TO EXCEL 
or KG 

c ****************************************************************** 

c 

DIMENSION XLAB2 ( 93 ) 

CHARACTER * 1 6 XLABl , XLAB2 
C 

DATA XLABl / * names 1 / 

DATA NSIMCH2 / 93 / 

C****************************************************************** 

c 

Cjcy INCLUDE ' dpguststats . inc ' 

c 

C****************************************************************** 

C 

DATA XLAB2 / 'names 1 , 

C 

c ******* 35 rms ******************************** 


★ 

1 DRMSTURBU ' , 

' DRMSTURBV ' , 

1 DRMSTURBW ' , 

* 

' DRMSTURBP • , 

' DRMSTURBQ ' , 

’ DRMSTURBR 1 , 

* 

1 DRMSTURBQ2 1 

, ' DRMSTURBR2 ' , 


* 

' DRMSFILU ' , 

1 DRMSFILV ' , 

1 DRMSFILW ' 7 

* 

' DRMSFILP ’ , 

' DRMSFILQ ' , 

' DRMSFILR 1 , 

★ 

' DRMSFILNU ' , 

' DRMSFILNV 1 , 

1 DRMSFILNW 1 , 

★ 

1 DRMSTURBUD ' , 

r 1 DRMSTURBVD 1 , 

’ DRMSTURBWD ' f 

* 

' DRMSFILUD ' , 

' DRMSFILVD ‘ , 

1 DRMSFILWD 1 , 

★ 

' DRMSFILUD2 ' , 

, 1 DRMSFILVD2 1 , 

' DRMSFILWD2 ' , 

* 

' DRMSTURALP ' , 

, ’ DRMSTURBBET ' 

, 

* 

1 DRMSTURALPD 1 

1 , ’ DRMSTURBBETD 

9 

★ 

n 

' DRMSFILA ’ , 

' DRMSFILB ' , 

1 DRMSFILAD ' , 

c ******* 37 means + 

r* 

2 rms +12 mil 

************** 

* 

' DMNTURBU ’ , 

' DMNTURBV ' , 

1 DMNTURBW 1 , 

* 

' DMNTURBP ’ , 

' DMNTURBQ ’ , 

’ DMNTURBR 1 , 

★ 

’ DMNTURBQ2 1 , 

' DMNTURBR2 ’ , 


* 

* DMNFILU 1 , 

' DMNFILV 1 , 

' DMNFILW ' # 

* 

’ DMNFILP 1 , 

' DMNFILQ ' , 

’ DMNFILR ' , 

* 

’ DMNFILNU 1 , 

1 DMNFILNV ' , 

' DMNFILNW 1 , 

* 

1 DMNTURBUD 1 , 

' DMNTURBVD ' , 

1 DMNTURBWD ' , 

* 

’ DMNFILUD 1 , 

' DMNFILVD 1 f 

' DMNFILWD 1 , 

* 

1 DMNFILUD2 ’ # 

' DMNFILVD2 ' , 

' DMNFILWD2 ' , 

* 

1 DMNTURBALP 1 , 

' DMNTURBBET ’ , 


* 

* DMNTURBALPD 1 

, ' DMNTURBBETD ' , 


★ 

1 DMNFILA 1 , 

‘ DMNFILB ' , 

1 DMNFILAD 1 , 

* 

’DRMSFILAD2 1 , 

' DRMSFILBD2 1 , 

' DMNFILAD 2 1 , 

★ 

' DRMSMILUK ’ , 

' DMNMILUK ' , 

' DRMSMILVK 1 , 

* 

' DMNMILVK ' , 

1 DRMSMILWK ’ , 

' DMNMILWK ' , 

* 

1 DRMSMILPK 1 , 

' DMNMILPK 1 , 

' DRMSMILQK ' , 

* 

1 DMNMILQK ' , 

* DRMSMILRK ' , 

1 DMNMILRK ' , 

* 

1 DRMSMILUDK ' , 

* DMNMILUDK 1 , 

' DRMSMILVDK ' f 

* 

' DMNMILVDK ' , 

' DRMSMILWDK ' , 

' DMNMILWDK ' / 


1 DRMSFILNP ' , 


1 DRMSFILBD ' , 


' DMNFILNP ' , 


1 DMNFILBD ' , 

' DMNFILBD2 1 , 
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C OPEN A DATA CHANNEL FOR WRITING THE GETDATA ASC FORMAT FILE. 

C 

OPEN (33, FILE= ' KGACSL . ASC ' , STATUS= ' UNKNOWN 1 ) 

C 

C WRITE OUT DATA IN TO FILE FOR USE IN GETDATA. 

C 

WRITE (33 format asc 2 . 1 ' 1 , / , ' ' nChans ' ' , t!4 , i3 ) ' ) NSIMCH2 -1 

WRITE (33 , ' (6al3) 1 ) XLAB1 , ( XLAB2 (I) , 1=2 , NSIMCH2 ) 

WRITE ( 33 , • ( ’ ' dataOOl ' ' ) 1 ) 

C 

WRITE ( 33 , 1 (6G13 .7) • ) T, 


* 

DRMSTURBU , 

DRMSTURBV, 

DRMSTURBW, 


* 

DRMSTURBP , 

DRMSTURBQ , 

DRMSTURBR , 


* 

DRMSTURBQ2 , 

DRMSTURBR2 , 



* 

DRMSFILU, 

DRMSFILV, 

DRMSFILW, 


* 

DRMSFILP, 

DRMSFILQ , 

DRMSFILR, 


* 

DRMSFILNU, 

DRMSFILNV, 

DRMSFILNW, 

DRMSFILNP t 

★ 

DRMSTURBUD, 

DRMSTURBVD , 

DRMSTURBWD , 


★ 

DRMSFILUD, 

DRMSFILVD, 

DRMSFILWD, 


★ 

DRMSFILUD2 , 

DRMSFILVD2 # 

DRMSFILWD2, 


★ 

DRMSTURALP, 

DRMSTURBBET, 



* 

DRMSTURALPD , 

DRMSTURBBETD , 



★ 

DRMSFILA, 

DRMSFILB , 

DRMSFILAD, 

DRMSFILBD , 

******* 37 means + 

2 rms + 12 mils ***********************’ 

* 

DMNTURBU , 

DMNTURBV , 

DMNTURBW , 


* 

DMNTURBP, 

DMNTURBQ, 

DMNTURBR , 


★ 

DMNTURBQ2 , 

DMNTURBR2 , 



* 

DMNFILU , 

DMNFILV, 

DMNFILW , 


* 

DMNFILP, 

DMNFILQ, 

DMNFILR, 


* 

DMNFILNU, 

DMNFILNV, 

DMNFILNW, 

DMNFILNP , 

* 

DMNTURBUD, 

DMNTURBVD , 

DMNTURBWD , 


* 

DMNFILUD, 

DMNFILVD, 

DMNFILWD , 


* 

DMNFILUD2 , 

DMNFILVD2 , 

DMNFILWD2 , 


* 

DMNTURBALP , 

DMNTURBBET , 



* 

DMNTURBALPD , 

DMNTURBBETD, 



* 

DMNFILA, 

DMNFILB , 

DMNFILAD, 

DMNFILBD , 

* 

DRMSFILAD2 , 

DRMSFILBD2 , 

DMNFILAD2 , 

DMNFILBD2 

★ 

DRMSMILUK, 

DMNMILUK, 

DRMSMILVK , 

DMNM I LVK , 

* 

DRMSMILWK, 

DMNMILWK, 

DRMSMILPK, 

DMNMILPK, 

* 

DRMSMILQK, 

DMNMILQK, 

DRMSMILRK, 

DMNMILRK , 

* 

DRMSMILUDK , 

DMNMILUDK , 

DRMSMILVDK, 

DMNMILVDK 

* 

DRMSMILWDK , 

DMNMILWDK 




CLOSE (33) 


RETURN 

END 

SUBROUTINE ASCFM (T, 


* 

TURBU, 

TURBV, 

TURBW, 

FILU , 

* 

FILV, 

FILW, 

TURBUD, 

TURBVD , 

* 

TURBWD , 

FILUD, 

FILVD , 

FILWD, 

* 

FILNU, 

FILNV, 

FILNW, 

FILNP, 

★ 

TURBP, 

TURBQ, 

TURBR , 


* 

FILP , 

FILQ, 

FILR , 

TURBALP 
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* 

* 

•k 


* 


TURBBET, 
FILB, 
TURBR2 , 
FILAD2, 
MILWK , 
MILUDK, 
NUMOUT , 


TURBALPD , 
F I LAD , 
FILUD2 , 
FILBD2 , 
MILPK , 
MILVDK , 
TSTP ) 


TURBBETD, 
FILBD , 
FILVD2 , 
MILUK, 
MILQK, 
MILWDK , 


FILA, 
TURBQ2 , 
FILWD2 , 
MILVK, 
MILRK , 


C 

q****************************************************************** 
C* SUBROUTINE ASCFM 

C* WRITES AN OUTPUT FILE IN THE GETDATA FORMAT ASC . 

c* ***************************************************************** 

c 

cjcy IMPLICIT REAL * 8 (a-h,o-z) 

c 

DIMENSION XLABEL2 ( 47 ) 

CHARACTER* 1 6 XLABELl , XLABEL2 
C 

cjcy was real *8 
C 

REAL MILUK, MILVK, MILWK, MILPK, MILQK, MILRK, 

* MILUDK, MILVDK, MILWDK 

C 

DATA XLABELl / 1 names 1 / 

DATA NSIMCH2 / 47/ 
c 

DATA XLABEL2 / 1 names 1 , 


* 

1 TURBU 1 , 

1 TURBV • , 

’ TURBW ’ , 

1 FILU 1 , 

* 

1 FILV 1 , 

* FILW 1 , 

1 TURBUD ’ , 

' TURBVD 1 

* 

' TURBWD 1 , 

1 FILUD 1 , 

’ FILVD’ , 

1 FILWD 1 , 

* 

' FILNU ' , 

1 FILNV 1 , 

* FILNW 1 , 

1 FILNP 1 , 

★ 

1 TURBP 1 , 

1 TURBQ 1 , 

’ TURBR 1 , 


* 

1 FILP ' , 

1 FILQ 1 , 

r FILR 1 , 

•TURBALP 

* 

' TURBBET 1 , 

’ TURBALPD ' , 

■ TURBBETD 1 , 

’FILA 1 , 

★ 

'FILB' , 

’ FILAD 1 , 

’ FILBD 1 , 

1 TURBQ2 1 

* 

1 TURBR2 1 , 

’ FILUD2 1 , 

’ FILVD2 1 , 

1 FILWD2 1 

* 

1 FILAD2 1 , 

’ FILBD2 ’ , 

1 MILUK 1 , 

’MILVK’ , 

* 

1 MILWK ’ , 

■MILPK 1 , 

'MILQK 1 , 

’MILRK 1 , 

* 

’MILUDK’ , 

1 MILVDK 1 , 

’MILWDK 1 / 



C 

C OPEN A DATA CHANNEL FOR WRITING THE GETDATA ASC FORMAT FILE. 

C 

IF ( (T .EQ. 0.0) ) THEN 

OPEN ( 3 , FILE= ' ACSL . ASC ' , STATUS= 1 UNKNOWN ' ) 

C 

C WRITE OUT DATA IN TO FILE FOR USE IN GETDATA. 

C 

WRITE (3, ' format asc 2 . 1 ' nChans tl4 , i3 )') NSIMCH2-1 

WRITE (3, ’ (6al3) ' ) XLABELl, (XLABEL2 (I) , 1=2 , NSIMCH2 ) 

WRITE (3, ' ( ' ' dataOOl ' ■ ) ■ ) 

INTERVA = NUMOUT 
END IF 
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c 


c 


IF (T .GE. 0.0 .AND. T . LT . TSTP) 
IF (INTERVA .EQ. NUMOUT) THEN 
WRITE ( 3 , * (6G13.7) 1 ) T, 


* 

TURBU , 

TURBV , 

TURBW, 

★ 

FILV, 

FILW, 

TURBUD , 

* 

TURBWD , 

FILUD, 

F I L VD , 

* 

FILNU , 

FILNV, 

FILNW, 

★ 

TURBP , 

TURBQ, 

TURBR , 

* 

FILP , 

FILQ, 

FILR, 

* 

TURBBET , 

TURBALPD, 

TURBBETD 

★ 

FILB , 

FILAD, 

FILBD, 

★ 

TURBR2 , 

FILUD2 , 

FILVD2, 

* 

FI LAD 2 , 

FILBD2 , 

MILUK, 

★ 

MILWK, 

MILPK , 

MILQK, 

* 

MILUDK, 

MILVDK, 

MILWDK 


INTERVA = 1 

ELSE 

INTERVA = INTERVA + 1 
END IF 
END IF 

IF ( T .GT. TSTP) CLOSE ( 3 ) 

RETURN 

END 


THEN 


FILU, 
TURBVD , 
FILWD, 
FILNP, 

TURBALP, 
FI LA, 
TURBQ2 , 
FILWD2 f 
MILVK , 
MILRK , 


29 



Gust amplitude, rad/sec 


Time History Plots 

Numerous 100-second runs were made with the GUSTMDL test program to produce turbulence sequences 
for analysis of the continuous, Tustin, and MIL STD turbulence models. Simulated aircraft velocities of 
100 ft/sec and 1000 ft/sec were used. Sample time history plots of the first 10 seconds of these sequences for 
each model are included in figures 1 through 3 for V = 100 ft/sec and in figures 4 through 6 for V = 1 000 
ft/sec. 



(a).- u-component 




(c).- w-component 





(e).- q-component (f).- r-component 

Figure 1. Sample time histories; V =100 ft/sec - continuous model. 


30 
















Gust amplitude, rad/sec 







Gust amplitude, rad/sec 










Gust amplitude, rad/sec 



Time, sec 


Time, sec 


(e).- q-component (f).- r-component 

Figure 5. Sample time histories; V =1000 ft/sec - Tustin discrete model. 
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(a).- u-component 







i r j ....... i x_i 

0 2 4 6 8 10 


Time, sec 


(e).- q-component (f) ~ r-component 

Figure 6. Sample time histories; V =1000 ft/sec - MIL STD discrete model. 


Power Spectral Densities 

Performance of the three implemented turbulence models was evaluated by computing the power spectral 
densities (PSD) of results of simulated turbulence from the GUSTMDL program. These PSD’s of the 
simulated turbulence components were then compared to the theoretical Dryden spectral model components. 
Equations used to compute the PSD for each model are listed in the section below, and the MATLAB m-files 
used to implement these equations and produce the PSD plots follow in the next section. 


Equations 

The two-sided theoretical power spectral densities for the continuous model were provided by the 
Government from the Dryden model (ref. 1 ) and are listed below as equations (36) through (41). 
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(J^T 1 

" r " / 

(36) 

s v {(o) - ^-0 v m - °* Tv 1 + 3(T v") 2 
2 2k (l + (V») 2 ) 

(37) 

■U°» =;*.(«»= *■ r - , 1+3(T - <a 2 \ i i 

2 2n (l + (t w (o) 2 ) 

(38) 

S p (a» = &£ 1 

1 + (v°) 

(39) 

V)- [( (<8/ t 

(40) 

*<«)- r , (<u/v)2 **<•) 
1 + [(3^M«/^)] 

(41) 


The theoretical power spectral densities for the Tustin model were provided by the Government and are 
listed below as equations (42) through (51). 


2n 


^ 2 ,, 

1+ COS| 

K)] 


* 1 

f 2 

+ C BL 

)+ 

( 2 
co u 

\ 

r*2 

C BL 

)cos(<yr v ) 




for i = v, w 


Kf of |a 2 + b 2 + c 2 + 2 (ab + bc)cos(a)T v ) + 2accos(2<or v )j 
d 2 +e 2 +f 2 + 2 (de + ef)cos((oT v ) + 2df cos(2coT v ) 


(43) 
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where 


®/ - ym 

K: = 


3 £0, 


In 

a = + <U ( - V 3 

b = 2(» ( - / a/3 
c = (0 -J a/3 - C flL 
d = (o>, + Cfl L ) 

« = 2(<»?-CSl) 

/ = ( w t -Ql) 


for / = v, vv 


(44) 




2 

c^M=^ Wv ) QwH 


where 


and 


where 


and 


p [ 

l+cos(ft>7 v ) 


* 1 

[oi 2 p + C BLp 

+ l 

( 2 r 2 \ 
\ ( °P- C BL P ) 

IcOS^^Ty) 


H, 


r^j 


2 2C 




i - cos(t»r v ) 


V 2 M 2 + + M^_ + 2M q+ Mq_ cos(coT v ) 


M q+ ~ 1 + T q C BL q 

M q - = l-T q C BLq ) 


G dr {(D)=H r [e ia>T ^G dv {(o) 


J_f ( iCOTy \ 

2 _ 2C BL r 

1 - cos(ft)T v ) 

H r \e V J 

V 2 

M 2 + + M 2 _ + 2 M r+ M r . 


M r + = 1 + X r C BL r 
M r _ = 1 - t r C BLr 


(45) 


(46) 


(47) 


(48) 


(49) 


(50) 


(51) 


The theoretical power spectral densities for the MIL STD model were provided by the Government and 
are listed below as equations (52) through (57). 
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where 


and 


where 


c <"(“)= frM^i 


KT; 


for / = u t v, w, p (52) 


1 + (l - a,-) 2 - 2(l - a, )cos(<»7^) 

H i [e j<oTv ) = < 7 ^'P- a i 

1 > l-d-aAe-W 


<J U y 2a,- 


for i = u, v, w, p (53) 


[l - (l — a/)cos(ft)7^,)] - 7(1 - a,)sin(<B7^) 


T v 

a u = — 

2T V 

a v= — 


2T V 

a w = — 


k W 


a P 


(54) 


#,(*) = 



for i,j = q,w or r,v 



2cos(fi)r v )] 

(55) 

(l + X 2 + y 2 ) - 2X(l + y)cos(<ar v ) + 2Ycos(2coT v ) 


5 («) 

for 1 , j = q,w or r, v 


Vj(z) 


k ikj(l-z~') 


(56) 

l ~( 2 ~ a i ~ a j)z~ l +(!-«,- 



*) 

1 - Xz~ l + Y Z ~ 2 
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and 



K = 


/r 




jfZv 

\ 


£- = (T 

lv w w 


4r„ 


V 

X = 2- cij - cij 

J' = (l- 0 ,)( , -a y )j 


for /,j = ?,wor r,v (57) 


Code 

The PSD equations listed in the previous section were coded into MATLAB m-file figl2dat.m. Time 
history sequences from the GUSTMDL simulation were input using the gdload utility avoiding conversion of 
input files prior to use of the MATLAB m-files. Execution of fig 1 2dat.m required assigning a value for the 
velocity V to correspond with the velocity used in the input file obtained from the GUSTMDL program. 
Using the gdwrite utility, the m-file produces output files in the asc2 format which is compatible with one of 
the plotting programs used by the DCB researchers. Source code for the m-file fig J 2dat.m is presented 
below. 


% Computes data for PSD plots f ile-f igl2dat .m 

% 

%********************** 

% First calc theoretical psd’s 
^********************** 

% 

% Snu = sampling function transform 

% Hxsq = magnitude in dB of transfer function for x-component 
% Sx = spectrum of q-, r-components 
% Tnu = sampling interval = 1/80 
% Lu = Lw = Lv = Lvw = Dryden scale length 
% sigma = sigu = sigw = sigv = std dev of turbulence 

Tnu = 0.0125; 

V = 100 . ; 
bwing = 37.42; 

Lu = 1750; 

Lwv = 1750; 

Lp = sqrt (Lwv*bwing) /2 . 6 ; 
s i gma - 5 . ; 
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sigp = 1 . 9/sqrt (Lwv*bwing) *sigma; 

Npts = 200; 

ft = logspace ( -2 , 2 , Npts ) ' ; 
w = 2 . *pi * f t ; 

tauu = Lu/V; 

tauwv = Lwv/V; 

taup = Lp/V; 

tauq = 4 . *bwing/ (pi*V) ; 

taur = 3 . *bwing/ (pi*V) ; 

Kwv = sigma A 2*tauwv/ (2*pi ) ; 

Ku = sigma^2 * tauu/pi ; 

Kp = sigp^2*taup/pi ; 

<£****★★***★***★★**★**** 

% Pre-allocate vectors 

<^* + ******************** 

Snu = zeros ( size (w) ) ; 

Swv = zeros (size (w) ) ; 

Su = zeros ( size (w) ) ; 

Sp = zeros ( size (w) ) ; 

Hqsq = zeros ( size (w) ) ; 

Hrsq = zeros ( size (w) ) ; 

Sq = zeros ( size (w) ) ; 

Sr = zeros ( size (w) ) ; 

Gdu = zeros (size (w) ) ; 

Gdwv = zeros ( size (w) ) ; 

Gdp = zeros ( size (w) ) ; 

Gdq = zeros { size (w) ) ; 

Gdr = zeros ( size (w) ) ; 

% 

Gmilu = zeros ( size (w) ) ; 

Gmilwv = zeros (size (w) ) ; 

Gmilp = zeros ( size (w) ) ; 

Gmilq = zeros ( size (w) ) ; 

Gmilr = zeros ( size (w) ) ; 

% Calculate psd's 

%********* + ************ 
for i = l:Npts, 

Snu(i) = 10 . *logl0 (sin (w(i ) *Tnu/2 . ) ^2/ (w(i) *Tnu/2 . ) ~2) ; 

N = 1. + 3 . * ( tauwv*w ( i ) ) ^2 ; 

D = (1 + ( tauwv*w ( i ) ) A 2 ) ^2 ; 

Swv(i) = 10 . *logl0 (Kwv*N/D) ; 

Su ( i ) = 10 . *logl0 (Ku/ ( 1 . + (tauu*w(i) ) * 2 ) ) ; 
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Sp ( i ) = 10 . *logl0 (Kp/ (1. + (taup*w(i) ) ~2) ) ; 


Hqsq(i) = 10 . *logl0 ( (w(i) /V) ~2/ (1 
Hrsq(i) = 10 . *logl0 ( (w(i) /V) "2/ (1 
Sq ( i ) = Hqsq ( i ) + Swv ( i ) ; 

Sr(i) = Hrsq(i) + Swv(i); 
end 
% 

^***************** 

% calc discrete theoretical psd * s 

<l*** + ************* 

% 

^***************** 

% u-component constants 
^***************** 

% 

wud = l./tauu; 

cfil = wud/tan (wud*Tnu/2 . ) ; 

sdnois = 1 . ; 

kud = sigma*sqrt ( tauu/pi ) ; 

cud = sdnois*sqrt (pi* (wud + cfil) 

ufacl = kud^2 *wud^2 ; 

udl = wud^2+cf i 1^2 ; 

ud2 = wud^2 -cf il^2 ; 

% 

% Mil constants - u 
au = Tnu/tauu; 
onemau = 1. - au; 
onemausq = onemau A 2; 

Tnusq = Tnu^2 ; 
sigusq=sigma^2 ; 

milunum= Tnusq/ (pi*tauu) *sigusq; 


% 

^***************** 
% v and w-component constants 
^***************** 

% 

wvdw = l./tauwv; 

% 

kwvd = sqrt (3 . / (2 . *tauwv*pi) ) ; 
wosr3 = wvdw/sqrt (3 . ) ; 
a = cfil + wosr3 ; 
b = 2 . *wosr3 ; 
c = wosr3- cfil; 
wnc = wvdw+cfil; 
d = wnc^2; 

e = 2 . * (wvdw A 2 -cf il^2 ) ; 
f = (wvdw-cf il ) ^2 ; 


( tauq*w ( i ) ) A 2 ) ) ; 
(taur*w ( i ) ) ^2 ) ) ; 


* ★ * ★ ★ 

★ ★ + ★* 

★ ★ + * * 

* * * + * 

) ; 
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wvfacl = kwvd A 2*sigma A 2 ; 
abcsum - a A 2+b A 2+c A 2; 
defsum = d A 2+e A 2+f A 2; 

% 

% Mil constants - w, v 
awv = 2 . *Tnu/ tauwv; 
onemawv =1. - awv; 
onemawvsq = onemawv A 2 ; 
milwvnum= Tnusq/ (pi*tauwv) *sigusq; 


% 

!£■*■*★★ + + ★★******★*★***** 
% p-component constants 
^***************** + + *** 
% 

wpd = l./taup; 

pcfil = wpd/ tan (wpd*Tnu/2 . ) ? 
mapt = (1. + pcfil*taup) ; 

mmpt = (1. - pcfil*taup) ; 
kpd = sigp*sqrt ( taup/pi ) ; 

% 

% 

% Mil constants - p 
ap = Tnu/taup; 
onemap = 1 . - ap ; 
onemapsq = onemap A 2 ; 
sigpsq = sigp A 2; 

milpnum= Tnusq/ (pi*taup) *sigpsq; 

%************★********* 
% q- component 

%********************** 

% 

wqd = l./tauq; 

cfilq = wqd/tan (Tnu/ (2 . *tauq) ) ; 
maqt = (1. + cfilq*tauq) ; 
mmqt = (1. - cfilq* tauq ) ; 
kqd = cfilq/V; 

% 

% 

% Mil constants - q 
aq = Tnu/ tauq; 
onemaq = 1 . - aq ; 

% onemawv defined above 
Xqw = 2 . - aq - awv; 

Yqw - onemaq* onemawv; 
milqwkl = 1. + Xqw A 2 + Yqw A 2; 
milqwk2 = 2.*Xqw*(l. + Yqw); 
milqwk3 = 2.* Yqw; 
kq = pi/ (4 . *bwing) ; 
kw = sigma*sqrt (2 . *awv) ; 

Tnuov2pi = Tnu/(2.*pi); 
kqsq=kq A 2 ; 
kwsq=kw A 2 ; 
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milqwnkl= Tnuov2pi*kqsq*kwsq; 

% milqwnkl= Tnuov2pi A 3 *kqsq*kwsq; 

% 

%********************* 
% r-component constants 
%********************* 


wrd = l./taur; 

cfilr = wrd/tan (Tnu/ (2 . *taur) ) ; 
mart = (1. + cfilr*taur) ? 
mmrt = (1. - cfilr*taur) ; 
krd = cfilr/V; 

% 

% Mil constants - r 
ar = Tnu/taur; 
onemar = 1. - ar; 

% onemawv, awv defined above 
Xrv =2. - ar - awv; 

Yrv = onemar* onemawv; 
milrvk4 = 1. + Xrv"2 + Yrv"2 ; 
milrvk5 = 2.*Xrv*(l. + Yrv); 
milrvk6 = 2.* Yrv; 
kr = pi/ (3 . *bwing) ; 

% kw defined above 
kv = kw ; 

% Tnuov2pi defined above 
krsq=kr^2 ; 
kvsq=kv^2 ; 

milrvnkl= Tnuov2pi*krsq*kvsq; 

% milrvnkl= Tnuov2pi^3 *krsq*kvsq ; 

% 

%********************* 
% Calculate psd’s 
% 

% u-component 

%********************* 

% 

for i = 1 : Npts , 

coswt = cos (w ( i ) *Tnu) ; 
cos2wt = cos ( 2 *w { i ) *Tnu ) ; 

UN = ufacl*(l. + coswt) ; 

UD = (udl + (ud2*coswt ) ) ; 

Gdu(i) = 10. *logl0 (UN/UD) ; 

% 

% Mil specs - u 
% 

milud = 1. + onemausq - 2 . *onemau*coswt ; 
Gmilu ( i ) = 10 . *logl0 (milunum/milud) ; 

% 

£********************* 

% v and w- component 

^********************* 



wvf ac2 
wvf ac3 
wvfac4 
wvf ac5 


2 . * (a*b + b*c)*coswt; 
2 . *a*c*cos2wt ; 

2 . * (d*e + e*f)*coswt; 
2 . *d* f *cos2wt ; 


WVN = wvfacl* (abcsum + wvfac2+ wvfac3); 

WVD = defsum + wvfac4 + wvfacS; 

Gdwv(i) = 10 . *loglO (WVN/WVD) ; 

% 

% 

% Mil specs - w,v 
% 

milwvd =1. + onemawvsq - 2 . *onemawv*coswt ; 
Gmilwv(i) = 10 . *logl0 (mi lwvnum/ milwvd) ; 

% 

%********************** 

% p-component 

%********************** 

% 

PN = kpd~2*2 . * (1 . + coswt ) ; 

PD = mapt ^2 + mmpt /v 2 + 2 . *mapt *irunpt * (coswt ) ; 
Gdp(i) = 10 . *logl0 (PN/PD) ; 

% 

% 

% Mil specs - p 
% 

milpd = 1. + onemapsq - 2 . *onemap*coswt ; 

Gmilp(i) = 10 . *logl0 (milpnum/milpd) ; 

% 

%**★******************* 

% q-component 

%********************** 

% 

QN = kqd^2 *2 . *( 1 . - coswt); 

QD = maqt A 2 + mmqt A 2 + 2 . *maqt *mmqt * ( coswt ) ; 
Gdq(i) = 10 . *logl0 (QN/QD) + Gdwv(i); 

% 

% 

% Mil specs - q 
% 

milqwnk2 = (2. - 2.*coswt); 
milqwnum = milqwnkl *milqwnk2 ; 

milqwd = milqwkl - milqwk2 *coswt + milqwk3 *cos2wt ; 
Gmilq(i) = 10 . *logl0 (milqwnum/milqwd) ; 

% 

%*★******************** 

% r-component 

%********************** 

% 
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RN = krd*2*2 . * (1. - coswt); 

RD = mart^2 + mmrt A 2 + 2 . *mart *mmrt* ( coswt ) ; 
Gdr(i) = 10 .* loglO (RN/RD) + Gdwv ( i ) ; 


% 

% Mil specs - r 


% 

% milqwnk2 defined above 
milrvnk2 = milqwnk2; 
milrvnum = milrvnkl *milrvnk2 ; 

milrvd = milrvk4 - mi IrvkS* coswt + milrvk6 *cos2wt ; 
Gmilr(i) = 10 . *logl0 (milrvnum/ milrvd) ; 

% 

end 

% 

%*******★********* + **** 

% Now calc MEASURED psd's of run gustrxxx . asc2 data 
^********************** 

% 

% First load data from run gustrxxx . asc2 
% 

gdload gustr32.asc2 
% 

% Now do continuous psd's 
% 

nfft = 1024; 
no v = 512; 
f s = 80; 

win = hanning ( 1024 ) ; 

SF = norm(win) ^2/sum(win) ^2; 

% 

[Puc,fm] = psd(turbu, nfft , fs, win, nov) ; 

Puc = 10 . *logl0 (Puc*SF) ; 

% 

[Pvc,fm] = psd (turbv, nfft , fs , win, nov) ; 

Pvc = 10 . *logl0 (Pvc*SF) ; 

% 

[Pwc,fm] = psd ( turbw, nfft , fs , win, nov) ; 

Pwc = 10 .* loglO ( Pwc*SF ) ; 

% 

[Ppc,fm] = psd ( turbp, nfft , fs , win, nov) ; 

Ppc = 10 .* loglO ( Ppc*SF) ; 

% 

[ Pqc , fm] = psd (turbq, nfft, fs, win, nov) ; 

Pqc = 10 .* loglO (Pqc*SF) ; 

% 

[Prc,fm] = psd(turbr, nfft, fs,win, nov) ; 

Prc = 10 . * loglO ( Prc*SF ) ; 

% 

% Now do discrete psd's 
% 

[ Pud, fm] = psd ( f ilu, nf f t , fs, win, nov) ; 
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Pud = 10 . *loglO ( Pud*SF) ; 

% 

[ Pvd, fm] = psd (filv,nfft, f s , win, nov) ; 

Pvd = 10 . * loglO ( Pvd*SF ) ; 

% 

[ Pwd, fm] = psd(filw,nfft , fs,win,nov) ; 

Pwd = 10.* log 1 0 ( Pwd* SF ) ; 

% 

[Ppd,fm] = psd(filp,nfft , fs, win, nov) ; 

Ppd = 10 . * loglO ( Ppd*SF ) ; 

% 

[Pqd # fm] = psd ( f ilq, nf f t , fs, win, nov) ; 

Pqd = 10 . *logl0 (Pqd*SF) ; 

% 

[Prd,fm] = psd ( filr,nf ft , fs, win, nov) ; 

Prd = 10 . *logl0 (Prd*SF) ; 

% 

% 

% Now do discrete mil spec psd ' s 
% 

[Pmilud, fm] = psd (miluk, nf f t , f s , win, nov) ; 

Pmilud = 10 . *logl0 (Pmilud*SF) ; 

% 

[Pmilvd,fm] = psd (milvk, nf ft , fs , win, nov) ; 

Pmilvd = 10 . * loglO ( Pmilvd*SF) ; 

% 

t Pmilwd, fm] = psd (milwk, nf f t , f s , win, nov) ; 

Pmilwd = 10 .* loglO ( Pmilwd*SF) ; 

% 

[Pmilpd, fm] = psd (milpk f nf f t , f s , win, nov) ; 

Pmilpd = 10 . *logl0 (Pmilpd*SF) ; 

% 

[Pmilqd, fm] = psd (milqk, nf f t , f s , win, nov) ; 

Pmilqd = 10 .* loglO ( Pmilqd*SF ) ; 

% 

[ Pmilrd, fm] = psd (milrk, nf f t , f s , win, nov) ; 

Pmilrd = 10 . *logl0 ( Pmilrd*SF) ; 

% Write results 

%******************* + ** 

% 

outvect = [ ' f t Snu Su Swv Hqsq Hrsq Sp Sq Sr Gdu Gdwv Gdp Gdq Gdr Gmilu 
Gmilwv Gmilp Gmilq Gmilr']; 

gdwrite { ' f igl2datt j . asc2 asc2 ' , outvect) 
gdwrite ( ’ f igl2datmj . asc2 asc2 ' , ' fm Puc,Pvc,Pwc 
Ppc, Pqc, Prc, Pud, Pvd, Pwd, Ppd, Pqd, Prd, Pmilud, Pmilvd, 

Pmilwd, Pmilpd, Pmilqd, Pmilrd 1 ) 

Plots 


PSD s of sample sequences of each turbulence component were produced using/zg/2dar./H for each of the 
three models. These data were input into a commercial plotting program, and the resulting plots are shown 
below in figures 7 through 12. These plots compare the sample PSD and theoretical PSD for each of the 
turbulence sequences that were shown in figures 1 through 6. 


46 







Power spectral density, dB P ower spectral density, dB Power spectral density, dB 



Frequency, Hz 



Frequency, Hz 


(a).- u-component 



Frequency, Hz 


(b).- v-component 



(c),- w-component 



Frequency, Hz 


(d).- p-component 



(e).- q-component (f).- r-component 

Figure 8. Power spectral densities, V = 100 ft/sec, Tustin model. 
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Power spectral density, dB Power spectral density, dB Power spectral density, dB 




Frequency, Hz 


(a).- u-component 


(b).- v-component 




(c).- w-component 


(d).- p-component 




(e).- q-component (0-* r-component 

Figure 9. Power spectral densities, V = 100 ft/sec, MIL STD model. 
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Power spectral density, dB Power spectral density, dB Power spectral density. dB 




(a).- u-component 


(b).- v-component 




(c).- w-component (d).- p-component 




(e).- q-component (f)*~ r-component 

Figure 11. Power spectral densities, V = 1000 ft/sec, Tustin model. 
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Power spectral density, dB Power spec tral density, dB p ower spectral density, dB 



Frequency, Hz 


Frequency, Hz 


(a).- u-component 



Frequency, Hz 

(c).- w -component 


(b).- v-component 



Frequency, Hz 

(d).- p-component 



Frequency, Hz 



Frequency, Hz 


(e).- q-component (f).- r-component 

Figure 12. Power spectral densities, V = 1000 ft/sec, MIL STD model. 
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Measured Statistics 


To evaluate the statistical accuracy of the simulated turbulence components extended runs were made with 
the GUSTMDL program. Each run was 1000 seconds long to reduce start-up effects and reduce the variance 
of sample statistics. A set of ten computer runs each using a different seed for the random noise generator 
was made using a velocity of 100 ft/sec and the parameter values identical to those for figures 1 through 6. A 
second set of runs using a velocity of 1000 ft/sec was generated. Sample statistics ( rms and mean values) 
were calculated for each component of turbulence for each of the implemented models by invoking the DRMS 
macro of the GUSTMDL program. 

Equations 

Sample statistics (mean and standard deviation ) were computed for each sequence for each 

computer run according to 




N 


k = 1 


i = w, v, w,/?, q, r; j = l,2,L,10 (58) 



Also computed were the mean Mg and standard deviation Eg of the sample means and sample 
standard deviations taken over the set of ten sequences for each component as follows: 



where g — fh^ or 6^ ; i = u,v,w,p,q,r. 

Results 

Results computed using equations (58) and (59) can be found in Tables 2 through 7. Sample means and 
standard deviations of each turbulence component for each of the ten 1000-second runs are shown for the 
continuous, Tustin, and MIL STD models. Also shown in the tables are the mean and standard deviation of 
the sample means and sample standard deviations taken over the set of ten sequences for each component 
using equations (60) and (61). 


53 



Table 2. Statistics for Ten 1000-Second Simulation Runs; 
Continuous Case; 1^=100 ft/sec 


Run 

Standard deviation 


u 

V 

w 

P 

Q 

r 

34 

4.95 

5.27 

4.90 

0.0371 

0.0212 

0.0242 

36 

4.80 

5.51 

5.83 

0.0375 

0.0214 

0.0245 

37 

5.27 

5.39 

4.82 

0.0373 

0.0209 

0.0239 

38 

4.84 

4.20 

5.16 

0.0366 

0.0209 

0.0236 

39 

5.27 

4.77 

4.51 

0.0379 

0.0214 

0.0243 

40 

5.24 

5.32 

4.76 

0.0394 

0.0207 

0.0236 

41 

5.11 

4.88 

5.00 

0.0366 

0.0210 

0.0245 

42 

5.46 

5.11 

4.78 

0.0372 

0.0209 

0.0240 

43 

5.36 

5.20 

5.07 

0.0368 

0.0208 

0.0236 

44 

4.40 

5.23 

5.24 

0.0360 

0.0215 

0.0243 

Mean 

5.07 

5.09 

5.01 

0.0372 

0.0211 

0.0241 

Std. Dev. 

0.321 

0.382 

0.358 

0.000912 

0.000273 

0.000369 

Run 

Mean 

- 



u 

V 

w 

P 

q 

r 

34 

0.296 , 

-0.112 

-0.463 

0.000106 

-6.42e-05 

6.54e-05 

36 

-1.23 

-0.951 

-0.388 

-0.000551 

3.45e-06 

2.08e-05 

37 

-0.261 

-0.740 

0.214 

-0.00264 

-2.90e-05 

-0.000109 

38 

0.259 

1.11 

-0.418 

-0.000202 

-5.56e-05 

5.06e-05 

39 

1.77 

-0.499 

-0.185 

-0.00267 

8.01 e-06 

-5.22e-06 

40 

Ml 

0.849 

0.606 

-6.02e-05 

2.72e-05 

1.67e-05 

41 

-0.419 

0.229 

-0.351 

0.000188 

-7.82e-05 

-3.97e-05 

42 

-1.09 

0.887 

-0.284 

-0.000544 

-4.67e-05 

-2.52e-05 

43 

0.425 

-1.26 

0.443 

-0.000526 

1.81e-05 

-6.10e-05 

44 

-0.660 

-0.258 

-0.923 

-0.00361 

-4. 1 5e-05 

5.71e-05 

Mean 

0.0204 

-0.0743 

-0.175 

-0.00105 

-2.58e-05 

-2.96e-06 

Std. Dev. 

0.950 

0.822 

0.463 

0.00138 

3.73e-05 

5.63e-05 
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Table 3. Statistics for Ten 1000-Second Simulation Runs; 
Tustin Case; V = 1 00 ft/sec 


Run 

Standard deviation 


u 

V 

w 

P 

q 

r 

34 

4.95 

5.23 

4.89 

0.0370 

0.0211 

0.0240 

36 

4.80 

5.44 

5.75 

0.0374 

0.0213 

0.0243 

37 

5.26 

5.35 

4.79 

0.0371 

0.0207 

0.0237 

38 

4.84 

4.18 

5.11 

0.0365 

0.0208 

0.0233 

39 

5.26 

4.74 

4.51 

0.0377 

0.0212 

0.0241 

40 

5.23 

5.27 

4.73 

0.0393 

0.0206 

0.0234 

41 

5.11 

4.84 

4.98 

0.0364 

0.0208 

0.0243 

42 

5.45 

5.08 

4.77 

0.0371 

0.0208 

0.0238 

43 

5.36 

5.15 

5.07 

0.0367 

0.0207 

0.0234 

44 

4.40 

5.21 

5.22 

0.0359 

0.0214 

0.0241 

Mean 

5.07 

5.05 

4.98 

0.0372 

0.0209 

0.0238 

Std. Dev. 

0.320 

0.373 

0.341 

0.000919 

0.000276 

0.000372 

Run 

Mean 


u 

V 

w 

P 

q 

r 

34 

0.295 

-0.0971 

-0.450 

0.000106 

-6.22e-05 

6.32e-05 

36 

-1.24 

-0.897 

-0.380 

-0.000551 

5.48e-06 

1.95e-05 

37 

-0.258 

-0.697 

0.205 

-0.00264 

-3.03e-05 

-0.000111 

38 

0.260 

1.06 

-0.408 

-0.000198 

-5.59e-05 

4.93e-05 

39 

1.77 

-0.474 

-0.176 

-0.00267 

7.45e-06 

-5.30e-06 

40 

1.11 

0.812 

0.585 

-6.05e-05 

2.74e-05 

1.45e-05 

41 

-0.413 

0.217 

-0.342 

0.000191 

-7.60e-05 

-3.90e-05 

42 

-1.08 

0.854 

-0.275 

-0.000548 

-4.70e-05 

-2.78e-05 

43 

0.422 

-1.18 

0.418 

-0.000532 

1.82e-05 

-6.23e-05 

44 

-0.661 

-0.239 

-0.877 

-0.00362 

-4.21e-05 

5.65e-05 

Mean 

0.0207 

-0.0646 

-0.170 

-0.00105 

-2.55e-05 

-4.22e-06 

Std. Dev. 

0.948 

0.781 

0.444 

0.00138 

3.70e-05 

5.62e-05 
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Run 


Standard deviation 



u 

V 

w 

P 

Q 

r 

34 

4.95 

5.29 

4.90 

0.0372 

0.0245 


36 

4.80 

5.47 

5.76 

0.0376 

0.0247 

1 

37 

5.27 

5.33 

4.85 

0.0374 

0.0241 


38 

4.84 

4.22 

5.11 

0.0368 

0.0241 


39 

5.26 

4.80 

4.59 

0.0380 

0.0246 

0.0282 

40 

5.24 

5.28 

4.81 

0.0395 

0.0239 

0.0274 

41 

5.11 

4.91 

4.93 

0.0367 

0.0242 

0.0284 

42 

5.46 

5.09 

4.80 

0.0373 

0.0241 

0.0278 

43 

5.36 

5.11 

5.03 

0.0369 

0.0240 

0.0274 

44 

4.40 

5.18 

5.23 

0.0361 

0.0248 

0.0282 

Mean 

5.07 

5.07 

5.00 

0.0373 

0.0243 

0.0279 

Std. Dev. 

0.320 

0.357 

0.320 

0.000914 

0.000312 

0.000421 

Run 


Mean 



u 

V 

w 

P 

q 

r 

34 

0.295 

-0.121 


0.000106 

r^6.13e-05 

5.09e-05 

36 

-1.24 

■ 

-0.368 

-0.000551 

5.42e-06 

3.08e-06 

37 

-0.258 

1 

0.203 

-0.00264 

-2.34e-05 

-0.000121 

38 

0.260 

1.10 

-0.423 

-0.000197 

-5.72e-05 

4,72e-05 

39 

1.77 

-0.491 

-0.191 

-0.00267 

9.76e-06 

-1.80e-05 

40 

1.11 

0.835 

0.617 

-6.16e-05 

1 .64e-05 

1.69e-05 

41 

-0.413 

0.232 

-0.341 

0.000191 

-6.41e-05 

-5.00e-05 

42 

-1.08 

0.866 

-0.295 

-0.000548 

-4.24e-05 

-2. 1 6e-05 

43 

0.422 

-1.27 

0.443 

-0.000532 

2.34e-05 

-6.00e-05 

44 

-0.661 

-0.262 

-0.930 

-0.00362 

-2.94e-05 

6.86e-05 


Mean 0.0206 -0.0819 -0.174 -0.00105 -2.23e-05 -8.37e-06 

Std-Dev. 0.948 0.821 0.464 0.00138 3.39e-05 5.83e-05 















Run 


Standard deviation 



u 

V 

w 

P 


r 

35 

5.07 

5.03 

5.07 

0.0373 

0.0207 

0.0238 

45 

5.17 

5.20 

5.22 

0.0372 

0.0207 

0.0238 

46 

5.16 

4.94 

4.98 

0.0373 

0.0206 

0.0237 

47 

5.10 

4.79 

5.05 

0.0370 

0.0207 

0.0237 

48 

5.28 

5.07 

5.16 

0.0373 

0.0208 

0.0240 

49 

5.31 

4.91 

4.92 

0.0376 

0.0208 

0.0237 

50 

5.26 

5.03 

5.03 

0.0370 

0.0206 

0.0240 

51 

5.22 

4.96 

5.03 

0.0372 

0.0207 

0.0239 

52 

4.92 

4.85 

5.01 

0.0370 

0.0206 

0.0238 

53 

4.91 

5.10 

5.20 

0.0368 

0.0207 

0.0238 

Mean 

5.14 

4.99 

5.07 

0.0372 

0.0207 

0.0238 

Std. Dev. 

0.141 

0.122 

0.0973 

0.000228 

8.38e-05 

0.000115 

Run 

Mean 


u 

V 

w 

P 

<7 

r 

35 

0.0393 


-0.159 

3.68e-05 

1.63e-08 

-4.04e-06 

45 

-0.417 


-0.115 

-0.000174 

-3.16e-06 

-2.90e-07 

46 

-0.0929 


0.0601 

-0.000824 

-8.60e-06 

-7.35e-06 

47 

0.0400 

0.362 

-0.147 

-5.01 e-05 

-9.49e-07 

-1.41e-06 

48 

0.551 

-0.157 

-0.0590 

-0.000864 

1.58e-06 

-9.09e-06 

49 

0.391 

0.271 

0.199 

-4.39e-05 

-6.69e-06 

-6.65e-06 

50 

-0.125 

0.0640 

-0.128 

5.42e-05 

4.89e-06 

-9.71e-06 

51 

-0.367 

0.270 

-0.104 

-0.000167 

-1.82e-07 

-3.87e-06 

52 

0.143 

-0.416 

0.143 

-0.000171 

8.47e-06 

-6.60e-06 

53 

-0.210 

-0.0693 

-0.303 

-0.00115 

2.91e-06 

8.59e-06 

Mean 

-0.00478 

-0.0266 

-0.0612 

-0.000335 

-1.72e-07 

-4.04e-06 

Std. Dev. 

0.308 

0.266 

0.152 

0.000437 

5.12e-06 

5.40e-06 















Run 


Standard deviation 



u 

V 

w 

P 

Q 

r 

35 

5.06 

5.02 

5.06 

0.0363 

0.0195 

0.0221 

45 

5.16 

5.18 

5.20 

0.0362 

0.0195 

0.0221 

46 

5.15 

4.93 

4.97 

0.0362 

0.0194 


47 

5.08 

4.77 

5.03 

0.0359 

0.0195 

0.0219 

48 

5.26 

5.05 

5.15 

0.0363 

0.0196 

0.0223 

49 

5.30 

4.90 

4.91 

0.0366 

0.0196 

0.0220 

50 

5.25 

5.02 

5.02 

0.0360 

0.0194 

0.0222 

51 

5.21 

4.95 

5.02 

0.0361 

0.0195 

0.0221 

52 

4.91 

4.84 

5.00 

0.0359 

0.0194 

0.0220 

53 

4.90 

5.08 

5.19 

0.0358 

0.0195 

0.0220 

Mean 

5.13 

4.97 

5.05 

0.0361 

0.0195 

0.0221 

Std. Dev. 

0.141 

0.122 

0.0977 

0.000235 

8.50e-05 

0.000120 

Run 

Mean 


u 

V 

w 

P 

<? 

r 

35 

0.0389 

-0.0219 


3.69e-05 

-5.64e-08 

-4. 1 1 e-06 

45 

-0.418 

-0.303 

EH 

-0.000174 

-3.61e-06 

-5.62e-08 

46 

-0.0919 

-0.266 

0.0605 

-0.000824 

-8.72e-06 

-7.21e-06 

47 

0.0404 

0.362 

-0.149 

-4.84e-05 

-7.64e-07 

-1.44e-06 

48 

0.551 

-0.157 

-0.0594 

-0.000865 

1.49e-06 

-9.28e-06 

49 

0.390 

0.271 

0.201 

-4.46e-05 

-6.83e-06 

-6.68e-06 

50 

-0.123 

0.0647 

-0.127 

5.50e-05 

4.66e-06 

-9.81 e-06 

51 

-0.365 

0.272 

-0.106 

-0.000168 

-1.10e-07 

-3.78e-06 

52 

0.142 

-0.416 

0.143 

-0.000173 

8.33e-06 

-6.40e-06 

53 

-0.211 

-0.0689 

-0.304 

-0.00115 

2.69e-06 

8.24e-06 

Mean 

-0.00473 

-0.0262 

-0.0615 

-0.000335 

-2.93e-07 

-4.05e-06 

Std. Dev. 

0.308 

0.267 

0.153 

0.000437 

5.12e-06 

5.34e-06 
















Run 


Standard deviation 



u 

V 

w 

P 

<7 

r 

35 

5.07 

5.06 

5.10 

0.0386 

0.0258 

0.0308 

45 

5.17 

5.19 

5.23 

0.0385 

0.0258 

: 0.0308 

46 

5.17 

4.97 

5.01 

0.0385 

0.0257 

0.0307 

47 

5.10 

4.81 

5.05 

0.0383 

0.0258 

0.0307 

48 

5.28 

5.06 

5.15 

0.0386 

0.0260 

0.0310 

49 

5.31 

4.93 

4.95 

0.0389 

0.0260 

0.0307 

50 

5.27 

5.06 

5.05 

0.0383 

0.0258 

0.0310 

51 

5.23 

4.98 

5.03 

0.0385 

0.0258 

0.0309 

52 

4.92 

4.88 

5.03 

0.0383 

0.0257 

0.0308 

53 

4.92 

5.10 

5.21 

0.0381 

0.0259 

0.0308 

Mean 

5.15 

5.00 

5.08 

0.0385 

0.0258 

0.0308 

Std. Dev. 

0.141 

0.112 

0.0906 

0.000227 

0.000102 

0.000132 

Run 

Mean 


u 

V 

vv 

P 


r 

35 

0.0387 

mzgm 

m 

3.71e-05 

-6.71e-07 

-3.69e-06 

45 

-0.418 

BH 

DU 

-0.000175 

-4.36e-06 

5.25e-07 

46 

-0.0919 

-0.264 

0.0600 

-0.000823 

-9.53e-06 

-4.69e-06 

47 


0.360 

-0.148 

-4.77e-05 

-2.1 le-07 

-1.92e-06 

48 


-0.157 

-0.0590 

-0.000865 

6.75e-07 

-9.44e-06 

49 

0.390 

0.270 

0.201 

-4.53e-05 

-7.01 e-06 

-8.05e-06 

50 

-0.123 

0.0644 

-0.126 

5.54e-05 

4.44e-06 

-9.36e-06 

51 

-0.365 

0.270 

-0.104 

-0.000167 

-3.17e-07 

-3.44e-06 

52 

0.142 

mm 

0.144 

-0.000173 

8.52e-06 

-6.00e-06 

53 

-0.211 

mm 

-0.303 

-0.00115 

2.65e-06 

7.13e-06 

Mean 

-0.00476 

-0.0263 

-0.0609 

-0.000335 

-5.82e-07 

-3.89e-06 

Std. Dev. 

0.308 

0.266 

0.153 

0.000437 

5.33e-06 

5.04e-06 



























Theoretical Statistics 


Theoretical, or expected, values for the sample statistics were desired for comparison with the measured 
values in Tables 2 through 7. This section presents the equations used for these calculations, the code 
implementing these equations, and results obtained in the form of tabulated data and plots. 

Equations 

The following equations defining the theoretical statistics were provided by the Government for 
implementation and computation. Equations (62) through (64) calculate the theoretical standard deviation of 
sample means of the various turbulence components. 
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Equations (65) through (68) define the theoretical standard deviation of sample standard deviations for 
each component. 
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where 
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%) = ~^-( R j- + R j+) + ^ L ( R i'~ + *‘> ) + ( 1/T . + 1/T “j ( R j- R, ~ + R J +R ‘ + ) 


2 2 

+ - 2 Rj-Rj-J - R] + + 2R j+ R j+t t) -?L( Rl 1 + R} + ) 

| 2 R i+ (R J+r T - R j+ )~ 2Rj_[Rj_ z T + > 4 (/? y+T /?, + + 

(l/T.+l/Ty ) 2 (l/T,.+l/T y ) 3 

3 ^ 4 

+ 7[v(W- 2 V) + *H(*H 7 ' +2,, ;-)]-f( ,i J-' t ^') 


( 68 ) 


In the above equations T , defined by T = ATT V , is the length of the turbulence sequence. (N is the 
number of samples in the sequence.) The variables CT, T, V, T v were defined previously in the section 
entitled Continuous Model. The variables /?,_ Rj_ R i+ Rj + Rj-r, Rj+ r are defined in the subsequent 
section discussing autocorrelation. 


Code 

Equations (62) through (68) were inplemented in the MATLAB m-file newA41.m. Comments were 
pHHpH to the File to identify equation numbers listed in the previous section. The variables Smhatuu, 
Smhatpp, Smhatqq , Smhatrr , and Smhatwv defined in file newA4J.m are the standard deviations of the 
sample means. The variables Ssighatuu, Ssighatpp, Ssighatqq , Ssighatrr, and Ssighatwv defined in file 
newA41.m are the standard deviations of sample standard deviations. 


%************************ file = newA41.m 

% Computes eq. A29, 62, 65 for u or p components 

% A37 -A43 , 64, 67, 62 for q or r components 

% A31,63, 66 for w or v components 

%************************************ 
% 

% Tnu = sampling interval = 1/80 

% Lu = Lw = Lv = Lvw = Dryden scale length 

% sigma = sigu = sigw = sigv = std dev of turbulence 

b = 37.42; 

Lu = 1750; 

Lwv = 1750; 

Lp = sqrt ( Lwv*b) /2 . 6 ; 
sigma = 5 . ; 

sigp = 1 . 9/sqrt (Lwv*b) *sigma; 

% add sigq uses sigw; sigr uses sigv; for 67 
%degs sigq = sqrt ( 0 . 05699*sigma A 2 ) ; 

%degs sigr = sqrt ( 0 . 07 667*sigma A 2 ) ; 
sigq = sqrt ( 1 . 73 6e-05*sigma A 2 ) ; 
sigr = sqrt ( 2 . 336e-05*sigma A 2 ) ; 

Tnu = 0.0125; 

Npts = 80000; 

% 
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% T=1000 sec, Npts=80000 

T=Npts*Tnu; 

% 

T2 = 2 . * T ; 

Tsg = T A 2 ; 

Tcube = T A 3; 

Tfour = T A 4 ; 

Ntau = 501; 

tau = linspace ( -10 , 10 , Ntau) ' ; 

tauu = Lu/V; 
tauwv = Lwv / V ; 
taup = Lp/V; 
tauq = 4 . *b/ (pi *V) ; 
taur = 3 . *b/ (pi*V) ; 

%********************** 

% Pre-allocate vectors (autocorrelation vectors) 
%********************** 

Rii = zeros ( si ze ( tau )) ; 

Rqq = zeros ( size ( tau) ) ; 

Rrr = zeros (size (tau) ) ; 

Ruu = zeros (size (tau) ) ; 

Rpp = zeros ( size (tau) ) ; 

Rwv = zeros (size(tau) ) ; 
%********************** 

% 

%********************** 

% 

% Logic to select u,w,p,q,r etc. for sigi, tauj , etc. 
% use variable selup, or selqr, etc for equations 
%*************★******** 

% 

%********************** 

% u, p Logic selection 
%*********★************ 

% 

%**★****★************** 

% Selected sigma and tau 
%********************** 

% 

if selup == ' u 1 
sigi = sigma; 
taui = tauu; 
elseif selup = = 'p* 
sigi = sigp; 
taui = taup; 

end 

% 

%********************** 

% eq. 69 eq. Rii autocorrelation function for u,p 
%********************** 

% 
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%********************** 

% tau - ’for loop' calculations 
%****★*★*************** 

% 

for n = 1 : Ntau 

% 

A29pwr = abs (tau(n) ) /taui; 

Rii(n) = sigi A 2*exp(-A29pwr) ; 

% 

if selup = = 'u' 

Ruu (n) = Rii (n) ; 
elseif selup == 'p' 

Rpp (n) = Rii (n) ; 

end 

end % end tau- loop for eq. 69 

% 

% 

%********************** 

% Compute expected values of sample statistics u or p components 
% 

%********************** 

% constants for 62 and 65 
^★********************* 

% eq. 62 

^********************** 

% for 62 and 65 
tisq = taui ^2 ; 

% constants for 62 
UPkl = (taui/T) ; 

UPk2 = ( tisq/Tsq) ; 

% 

^★********************* 

% eq. 62 for sigma (mean-hat ) 
%*******★************** 

% 

Smhat46 = sigi*sqrt (2 . * (UPkl - UPk2 ) ) ; 

% 

%***********★********** 

% eq. 65 

^********************** 

% constants for 65 and use tisq defined above 
ticube = taui ^3; 
tifour = taui^4; 

% 

UPk4=taui/T2; 

UPk5=9 . *tisq/ (4 . *Tsq) ; 

UPk6=4 . *ticube/ (Tcube ) ; 

UPk7 - 2 . *tif our/ (Tfour ) ; 

% 

^********************** 

% eq. 65 for sigma ( sigma-hat ) 

% 

Ssighat53 = sigi *sqrt (UPk4 - UPk5 + UPk6 - UPk7 ) ; 

% 
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%********************** 

% 

%********************** 

% Set all u or p calculations 
%*********★★*********** 

% 

if selup = = 'u' 

% Ruu (n) = Rii (n) ; 

Smhatuu = Smhat46; 

Ssighatuu = Ssighat53; 

Smhatpp = ’undef'; 

Ssighatpp = 'undef'; 
elseif selup == 'p' 

% Rpp (n) = Rii (n) ; 

Smhatpp = Smhat46; 

Ssighatpp = Ssighat53; 

Smhatuu = 'undef'; 

Ssighatuu = 'undef' ; 

end 

% 

%**** + + + ★*★**★** + * + **** 

% end u, p Logic selection and calculations 
%********★**********★*★ 

% 

% 

% 

%********************** 

% Start q, r Logic selection 
%********************** 

% 

% Selected sigma and tau i,j = q,w or r,v 
%*******★************** 

% 

if selqr == 'q' 
sigi = sigq; 
sigj = sigma; 
taui = tauq; 
tauj = tauwv; 
elseif selqr == 'r' 
sigi = sigr; 
sigj = sigma; 
taui = taur; 
tauj = tauwv; 

end 

% 

%********************** 

% 

% eq. 71 - 77 constants for Rii autocorrelation function 
% for q,r sigma (mean-hat ) 

% 

%*★*** ***************** 

% for eq. 64 also- use tisq, tjsq 
tisq = taui ^2 ; 
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tjsq = tauj A 2; 

% 

tauij = tauj*taui; 

% 

oneovtaui - l./taui; 
oneovtauj = l./tauj; 

RAk3 = oneovtaui+oneovtauj ; 

RAkl2 = oneovtaui /RAk3 ; 

RAk4 = 1 . / (2 . * tauij ) ; 

RAk3sq =RAk3^2; 

RA4ov3sq = RAk4/RAk3sq; 

RAk5 = oneovtaui -oneovtauj ; 

RAkl3 = oneovtaui /RAk5 ; 

RAkSsq =RAk5 A 2; 

RA4ov5sq = RAk4 /RAkSsq; 

oneovtauisq = l./tisq; 
oneovtaujsq = l./tjsq; 

RAk6 = oneovtaui sq- oneovtauj sq; 

RAkll = oneovtauisq/RAk6 ; 

RAk7 = RAk4*oneovtaui; 

RAk5by3sq = RAk5*RAk3sq; 

RA5sqby3 = RAk5sq*RAk3; 

RAk8 = 1 . / (2 . *tauj ) ; 

RAk4ov3 = RAk4/RAk3; 

RAk4ov5 = RAk4/RAk5; 

RAk7ov6 = RAk7/RAk6; 

RAk9 = 1 . / (2 . *taui ) ; 

RAk9ov5 = RAk9/RAk5; 

RAk9ov3 = RAk9/RAk3; 

RAklO = 1 . / (4 . *tauij ) ; 

% 

%RAk3cube used in Sk7den eq. 68 
RAk3cube =RAk3~3; 

% not - used RA3by5 = RAk3*RAk5; 

% 

^********* 

% Eq 72 

^********* 

Rjneg = 1. - RAkl2 + RA4ov3sq - RAkl3 - RA4ov5sq 

+ RAkll - RAk7 /RAk5by3sq + RAk7 /RA5sqby3 ; 

%***★★*★** 

% Eq 73 

%**★*★★*** 

Rjnegt = RAk8 - RAk4ov3 - RAk4ov5 + RAk7ov6; 
%********* 
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% Eq 74 

Rineg = RAk9ov5 -RAk9ov3 + RAklO /RAk5 sq + RAklO/RAk3sq; 

^★★******-Ar 

% Eq 75 

<£**★*•*★*** 

Rjpos =1. - RAkl3 - RA4ov5sq - RAkl2 + RA4ov3sq . . . 

+ RAkll + RAk7 /RASsqby 3 - RAk7/RAk5by3sq ; 

% Eq 7 6 

%***★**★★* 

Rjpost = - RAk8 + RAk4ov5 + RAk4ov3 - RAk7ov6; 

%********* 

% Eq 77 

Ripos = - RAk9ov3 + RAk9ov5 + RAklO/RAk3sq + RAklO /RAkSsq; 
% 

Ripsq = Ripos /v 2 ; 

Rinsq = Rineg /V 2; 

Rjpsq = Rjpos~2; 

Rjnsq = Rjneg^2; 

Rjptsq = Rjpost ^2 ; 

Rjntsq = Rjnegt /V 2; 

% 

%********************** 

% eq. 71 eq. Rii autocorrelation function for q,r 
%******************★*** 

% 

%********★★************ 

% start tau - 'for loop' for 71 
%**********★***★******* 

% 

% constant for 71 

Riikl = sig j ^2 / ( V^2 *tisq) ; 

% 

% 

for n = l:Ntau 


% 

% 


% 


% 


A37pwrl = tau(n)/tauj; 
A37pwr2 = tau(n)/taui; 


tau < 0 

if tau(n) 
Rii (n) 


else 

tau(n) >= 0 

Rii (n) 


end 


< 0 

= Riikl * (Rjneg*exp (A37pwrl ) ... 

+ Rjnegt*tau (n) *exp ( A37pwrl ) + Rineg*exp { A37pwr2 ) ) ; 


= Riikl* (Rjpos*exp ( -A37pwrl ) ... 

+ Rjpost * tau (n) *exp ( -A37pwrl ) + Ripos*exp ( -A37pwr2 ) ) 
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if selqr = = 'q' 

Rqq (n) = Rii (n) ; 
elseif selqr == 'r' 

Rrr (n) = Rii ( n) ; 

end 

% 

end % end tau-loop for eq . 71 

% 

%********★************* 

% 

% Compute expected values of sample statistics for q or r components 
% 

%********************** 

% eq . 64 

^★★★★★-fr*********"******* 

% 

% constants needed for i,j = q,w or r , v 
% for 64 and tisq, tjsq defined above 
Tti = T*taui; 

Tt j = T*tauj ; 
tjcube = tauj A 3; 

% for 68 also use tjcube, Tti, Ttj defined above 
tjfour = tau j A 4 ; 

% 

%ticube = taui A 3; 

%tifour = taui A 4; 

£********************** 

% 64 use tisq # tjsq calcs defined above 

% constants for 64 
% 

QRkl = sigj / (V*Tti) ; 

QRk2 = Ttj*(Rjneg + Rjpos) ; 

QRk3 = Tti*(Rineg + Ripos) ; 

QRk4 = T*t jsq* (Rjpost - Rjnegt); 

QRk5 = tjsq* (Rjpos + Rjneg) ; 

QRk6 = tisq* (Ripos + Rineg) ? 

QRk7 = 2 . *tjcube* (Rjnegt - Rjpost); 

QRk2toQRk7 = (QRk2 + QRk3 + QRk4 - QRk5 - QRk6 + QRk7 ) ; 

% 

%★************★******** 

% eq. 64 for q,r sigma (mean-hat ) 

% 

Smhat 50 = QRkl*sqrt (QRk2toQRk7) ; 

% 

^★********************* 

% 

%★********************* 

% eq. 68 constants for S00 
% 

Ski = Ttj/2 . * (Rjnsq + Rjpsq) ; 

Sk2 = Tti/2 . * (Rinsq + Ripsq) ; 

Sk3 = T2/RAk3* (Rjneg*Rineg + Rjpos*Ripos) ; 

Sk4 = t j sq/4 .*( -Rjnsq - 2 . *Rjneg*Rjnegt*T - Rjpsq + 2 . *Rjpos*Rjpost*T) 
Sk5 = tisq/ 4 .* (Rinsq + Ripsq); 



Sk6num = 2 . *Ripos* (Rjpost *T - Rjpos) - 2 . *Rineg* (Rjnegt*T 
Sk6den = RAk3sq; 

Sk7num = 4 . * (Rjpost*Ripos + Rjnegt*Rineg) ; 

Sk7den = RAk3cube; 

Sk8 = (Rjpost* (Rjpost*T - 2 . *R jpos ) ) ; 

Sk9 = (Rjnegt* (Rjnegt*T + 2.*Rjneg)); 

SklO = t jcube/4 . * (Sk8 + Sk9); 

Skll = 3 . *t j four/ 8 . * (Rjntsq + Rjptsq) ; 

% 

%**************, If******* 

% eq. 68 calculation - S00 constant 

% 

S00 = Ski + Sk2 + Sk3 + Sk4 - Sk5 + Sk6num/Sk6den . . . 

+ Sk7num/Sk7den + SklO - Skll ; 

% 

%*************** + + + + + + + 

% 

%************.^ ri(f + 1 n r5itr + ^ i ^^^ 

% eq. 67 
% 

% constants for 67 

QRk8 = sig j ^2 / (sigi* (V*taui) ^2*T*sqrt (2. ) ) • 

% 

.fc*****************^*^*^ 

% eq. 67 for sigma ( s igma-hat ) 

% 

Ssighat 62 = QRk8*sqrt (S00 - 2 . * ( QRk2toQRk7 ) A 2 ) ; 

% 

%*********************^ 

% 

%********************** 

% Set all q or r calculations 

%************ + **4' i ' il . i '1'i l .' k 

% 

if selqr == 'q' 

% Rqq (n) = Rii (n) ; 

Smhatqq = Smhat50; 

Ssighatqq = Ssighat62; 

Smhatrr = 'undef'; 

Ssighatrr = 'undef'; 
elseif selqr == ’r' 

% Rrr (n) = Rii (n) ; 

Smhatrr = SmhatSO; 

Ssighatrr = Ssighat62; 

Smhatqq = 'undef'; 

Ssighatqq = 'undef'; 

end 

% 

%************* ir + i ' ic + i ' :ki ' i ' 

% w or v calculation 

%********************** 

% 

%********************* + 


Rjneg) 
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% eq.70 Rwv autocorrelation function for v,w 
^************* + ******** 

% 

^********************** 

% tau - 'for loop' calculations 

C^********************** 

% 

for n = l:Ntau 

% 

abstau = abs(tau(n)); 

A31pwr = abstau/tauwv; 

WVkl = 3 . *sigma / '2 /tauwv; 

WVk2 = tauwv/ 3.; 

WVk3 = ( 1 . /6 . ) *abstau; 

Rwv (n) = WVkl* (WVk2 *exp ( -A3 Ipwr ) - WVk3 *exp ( -A3 lpwr ) ) ; 

% 

end % end tau-loop for eq. 70 

% 

^********************** 

% eq. 63 for v, w sigma (mean-hat ) 

% 

Smhat48 = sigma*sqrt ( tauwv/T) ; 

% 

^********************** 

% eq. 66 for v, w sigma ( s igma-hat ) 
^********************** 

% 

WVk4 = 13 . *tauwv/ ( 4 . *T) ; 

WVk5 = 83 *tauwv /N 2 / ( 8 . *Tsq) ; 

Ssighat55 = ( sigma/2 .) *sqrt (WVk4 - WVk5); 

% 

Smhatwv = Smhat 4 8 ; 

Ssighatwv = Ssighat55; 

% 

% Display expected values of sample statistics 
% 

Smhatuu 

Smhatpp 

Smhat qq 

Smhat rr 

Ssighatuu 

Ssighatpp 

Ssighatqq 

Ssighatrr 

Smhatwv 

Ssighatwv 


Autocorrelation 

The theoretical autocorrelation functions were needed to compute some of the theoretical statistics 
discussed previously. The functions were also used for comparison with correlation functions computed from 
the sample sequences produced with the GUSTMDL program. 
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Equations 


The autocorrelation function (t), i = u, p, for the u- and p-components can be found from the 
Fourier transform of the spectral density as follows: 
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The autocorrelation function for the w-component is 
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Since the spectra of the v-component and w-component are the same, the autocorrelation function /? vy (r) 
is also given by equation (70). For the q- and r-components 
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These equations were coded in the MATLAB m-file newA47.m presented in the previous section. 
Comments were added to the code to refer to equation numbers listed in this section. The variables Ruu, 
Rpp, Rqq , Rrr , and /?wv are the autocorrelation functions from equations (69) through (71) for the turbulence 
components. The variables Rjneg , Rjnegt, Rineg , Rjpos, Rjpost , and /?//7os in file newA4J.m are the 
constants defined in equations (72) through (77). 


Code in MATLAB file A41plts.m was generated to plot the theoretical autocorrelation functions 
discussed above that are calculated in file newA4J.m. The MATLAB files newA4Lm and A4iplts.m should 
be executed sequentially to be able to plot theoretical autocorrelation functions. The file A41plts.m also 
computes and plots the sample autocorrelation functions of the turbulence sequences produced by the 
GUSTMDL program. These autocorrelation functions were computed using the MATLAB function xcorr 
from the Signal Processing Toolbox. The code for A41plts.m follows. 


% ************ ****** ***** file A41plts.m 

% 

% autocorrelation plots 

% 

whitebg 

% 

% Theorectical autocorrelation plots calculated in newA41.m 

% 

% Select u or p component 
% 

if selup == ' u' 

% 

plot ( tau, Ruu) 

title (' Autocorrelation Ruu V100’) 
xlabel ( ' tau, sec ' ) 
ylabel ( 'Autocorrelation ' ) 

% 

%pause 
%print Ruu 
% 


71 


elseif selup == 'p' 

% 

plot ( tau ; Rpp) 

title ( 'Autocorrelation Rpp V100’) 
xlabel ( ' tau, sec ' ) 
ylabel ( 'Autocorrelation 1 ) 

% 

%pause 
%print Rpp 
% 

end 

% 

% Select q or r component 
% 

if selqr == 'q' 

% 

plot ( tau , Rqq) 

title (’ Theoretical autocorrelation Rqq V100 1 ) 
xlabel ( 1 tau, sec 1 ) 
ylabel ( 'Autocorrelation' ) 

% 

%pause 
%print Rqq 
% 

elseif selqr == 'r 1 

% 

plot ( tau , Rrr ) 

title ( ' Theoretical autocorrelation Rrr V100 1 ) 
xlabel ('tau, sec') 
ylabel ( 'Autocorrelation' ) 

% 

%pause 
%print Rrr 
% 

end 

% 

% Select w or v component 
% 

plot (tau, Rwv) 

title (' Theoretical autocorrelation Rwv V100 ' ) 
xlabel ( ' tau, sec ' ) 
ylabel ( 'Autocorrelation' ) 

% 

%pause 
%print Rwv 
% 

%******** ******** ******** ******** 


Experimental autocorrelation plots from GUSTMDL 



if V==100. 

% 


gdload ( ' gustr54sigs . cmp3 ' ) 

% 
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nx = f ix( length (tx) /2) + 1; 
tplt = tx (nx-800 :nx+800 ) ; 
Rwcplt = Rvvc (nx- 800 : nx+800 ) 
RwTplt = RwT (nx-800 : nx+800 ) 
RvvMplt = RvvM (nx- 800 : nx+800 ) 

% Plot w or v component 


% 

plot (tplt, Rwcplt, 'k-' , tplt, RwTplt, 'y-' , tplt, RvvMplt, 'm-' , . . . 
tau, Rwv, ' k- - ' ) 

% 

xlabel(’tau, sec') 
ylabel ( 'Autocorrelation' ) 
grid 

legend ( ' Rvvc 1 , 1 RwT 1 , 1 RwM ' , ' Rwv ' ) 

% 

if V==l 00 . 

% 

t it le ( 1 Experiment al vs. Theoretical autocorrelation Rw V100') 
axis ( [ -2 . 2 . -5 . 25 . ] ) 

% 

elseif V==1000. 

% 

title (' Experimental vs. Theoretical autocorrelation Rw V1000') 
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axis ( [-2. 2. -10 . 25. ] ) 

end 

elseif selvw=- ' w 1 

[Rwwc] = xcorr (turbw, turbw, ’biased'); 

[ RwwT ] = xcorr (turbw, turbw, 'biased'); 

[RwwM] = xcorr ( turbw, turbw, 'biased'); 

tx = ( [1: length (Rwwc) ] 1 - f ix ( length (Rwwc ) /2 ) - l)*dt; 

% 

% Select middle +-10 seconds 
% 

nx = fix ( length ( tx) / 2 ) + 1; 

tplt = tx (nx- 800 : nx+800 ) ; 

Rwwcplt = Rwwc (nx-800 : nx+ 800 ) ; 

RwwTplt = RwwT (nx-800 :nx+800) ; 

RwwMplt = RwwM (nx-800 :nx+800) ; 

% Plot w or v component 
% 

% 

plot (tplt , Rwwcplt , ' k- ' , tplt , RwwTplt , 'y- ' , tplt , RwwMplt , ' m- ' , . . . 
tau , Rwv, ' k- - ' ) 


% 

xlabel ( 1 tau, sec ’ ) 
ylabel ( 'Autocorrelation' ) 
grid 

legend ( ' Rwwx ' , ' Rwv ' ) 

% 

if V==100 . 

% 

title ( 1 Experimental vs. Theoretical autocorrelation Rww V100') 
axis ( [-2. 2. -5 . 25. ] ) 

% 

elseif V= = 1 0 0 0 . 

% 

title (' Experimental vs. Theoretical autocorrelation Rww V1000') 
axis ([-2. 2. -10. 25.]) 

end 

% 

pause 

end 

% 

% 

%****★*★* ******** ******** ******** 

% Select u or p component 

%******** ******** ******** ******** 

% 

if selup -= ’u' 

% 

[Ruuc ] = xcorr ( turbu , turbu, 'biased'); 

[ RuuT ] = xcorr ( turbu , turbu, 'biased'); 

[RuuM] = xcorr (turbu, turbu, 'biased'); 
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Ruucplt = Ruuc (nx-800 : nx+800 ) 
RuuTplt = RuuT (nx-800 : nx+800 ) 
RuuMplt = RuuM (nx-800 : nx+800 ) 


% 

plot ( tplt , Ruucplt , 1 k- ' , tplt , RuuTplt , 'y- 1 , tplt , RuuMplt , 1 m- 1 , . . . 
tau, Ruu, 1 k-- 1 ) 

% 

xlabel(’tau, sec*) 
ylabel ( 'Autocorrelation 1 ) 
grid 

legend ( ' Ruuc ' , ' RuuT 1 , ' RuuM ' , ' Ruu 1 ) 

% 

if V==100. 

% 

title ( 'Experimental vs. Theoretical autocorrelation Ruu V100') 
axis ([-2. 2. -5. 35.]) 

% 

elseif V==1000. 

% 

title (' Experimental vs. Theoretical autocorrelation Ruu V1000') 
axis ([-2. 2. -10. 30.]) 

end 

% 

pause 

% 

elseif selup == 1 p * 

% 

[ Rppc ] = xcorr(turbp, turbp, ’biased'); 

[ RppT ] = xcorr ( turbp , turbp, ’biased’); 

[RppM] = xcorr (turbp, turbp, ’biased'); 

Rppcplt = Rppc (nx-800 : nx+800 ) ; 

RppTplt = RppT (nx-800 :nx+800 ) ; 

RppMplt = RppM (nx-800 : nx+800 ) ; 

% 

plot (tplt , Rppcplt , ’ k- ’ , tplt , RppTplt, ’y~ 1 , tplt , RppMplt , ’ m- ’ , . . . 
tau , Rpp , 1 k- - ‘ ) 


xlabel ( 1 tau , sec 1 ) 
ylabel ( ’Autocorrelation’ ) 
grid 

legend ( ’ Rppx ’ , ’ Rpp ’ ) 

% 

if V==100. 

% 

title ( 1 Experimental vs. Theoretical autocorrelation Rpp V100’) 
axis ([-2. 2. -2.e-4 14.e-4]) 

% 

elseif V= = 1 00 0 . 
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title (' Experimental vs. Theoretical autocorrelation Rpp V1000') 
axis ([-2. 2. -2 . e-4 14.e-4]) 

end 


% 

pause 

% 

end 

% 

%******** ******** ******** ******** 

% 

% Select q or r component 

%******** ******** ******** ******** 

% 

if selqr == ’ q 1 

% 

[ Rqqc ] = xcorr(turbq, turbq, 'biased'); 

[ RqqT ] = xcorr (turbq, turbq, 'biased'); 

[RqqM] = xcorr(turbq, turbq, 'biased'); 

Rqqcplt = Rqqc (nx-800 :nx+800) ; 

RqqTplt = RqqT (nx-800 : nx+800 ) ; 

RqqMplt = RqqM (nx-800 : nx+800 ) ; 

% 

plot (tplt, Rqqcplt , ’ k- 1 , tplt, RqqTplt, ’y- • , tplt , RqqMplt , 'm- ' , . . . 
tau, Rqq, ’ k-- 1 ) 

% 

xlabel ( ’ tau, sec 1 ) 
ylabel ( 'Autocorrelation 1 ) 
grid 

legend ( ' Rqqc ’ , 1 RqqT ' , ' RqqM ' , ' Rqq ' ) 

% 

if V==100 . 

% 

title (' Experimental vs. Theoretical autocorrelation Rqq V100') 
axis ([-2. 2. -l.e-4 5. e-4]) 

% 

elseif V==1000. 

% 

title (' Experimental vs. Theoretical autocorrelation Rqq V1000') 
axis ([-2. 2. -l.e-4 5. e-4]) 

end 

% 

% 

pause 

% 

elseif selqr == 'r' 

% 

[Rrrc ] = xcorr (turbr, turbr, 'biased'); 

[RrrT] = xcorr(turbr, turbr, ’biased’); 

[RrrM] = xcorr (turbr, turbr, 'biased'); 
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Rrrcplt = Rrrc (nx-800 : nx+800 ) ; 
RrrTplt - RrrT (nx-800 : nx+800 ) ; 
RrrMplt = RrrM (nx-800 :nx+800) ; 


% 

plot (tplt, Rrrcplt, • k- ' f tplt , RrrTplt , 'y- 1 , tplt # RrrMplt , 'm- 1 , . - . 
tau, Rrr , 1 k-- 1 ) 

% 

xlabel('tau, sec') 
ylabel ( 'Autocorrelation' ) 
grid 

legend ( 1 Rrrc ' , 1 RrrT ' , ’ RrrM ' , ' Rrr ' ) 

% 

if V==100 . 

% 

title (’ Experimental vs. Theoretical autocorrelation Rrr V100') 
axis ([-2. 2. -l.e-4 6.e-4]) 

% 

elseif V==1000. 

% 

title (' Experimental vs. Theoretical autocorrelation Rrr V1000') 
axis ([-2. 2. -l.e-4 6.e-4]) 

end 


% 

pause 

% 

end 

if selup == 'u' & selvw == 'V & selqr == 'q' 
outvectl = ['tau Ruu Rwv Rqq tplt Ruucplt RuuTplt RuuMplt ' ]; 

outvect2 = [ 1 Rvvcplt RvvTplt RwMplt Rqqcplt RqqTplt RqqMplt ' ]; 

outvect = [outvectl outvect2] ; 
gdwrite ( 1 cor55uvq . asc2 asc2 1 , outvect ) 

elseif selup = = ’p’ & selvw == 'w' & selqr == 'r* 
outvectl = ['tau Rpp Rwv Rrr tplt Rppcplt RppTplt RppMplt 1 ]; 

outvect2 = [ ' Rwwcpl t RwwTplt RwwMplt Rrrcplt RrrTplt RrrMplt' ]; 

outvect = [outvectl outvect2] ; 
gdwrite ( ' cor55pwr . asc2 asc2 ' , outvect ) 

end 


Plots 

Example plots of autocorrelation functions produced by A4Jplts.m for V= 100 ft/sec and 
V= 1000 ft/sec are found in this section. The sample autocorrelation functions are for 819-sec turbulence 
sequences produced using GUSTMDL. 
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Autocorrelation Autocorrelation 



Lag, see 



Lag, sec 


(a).- u-component 


(b).- v-component 




Lag, sec 


(c).- w-component 


(d).- p-component 




Lag, sec Lag, sec 

(e).- q-component (f).- r-component 

Figure 13. Autocorrelation functions for V = 100 ft/sec. 
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Autocorrelation 








Aircraft Simulation 


The turbulence models and implementations discussed in previous sections of this report were developed 
for use in aircraft simulations, primarily nonlinear six-degree-of-freedom simulations. The models were 
installed in a nonlinear six-degree-of-freedom simulation of the High Alpha Research Vehicle (refs. 3 and 4) 
for evaluation and comparison among the three models. The HARV simulation is written in ACSL with 
FORTRAN subroutines to implement input/output, sensor models, and control laws. The equations of 
motion are integrated in the ACSL DERIVATIVE BLOCK (ref. 5). 

Code 

Much of the code installed in the HARV simulation to implement the turbulence models was taken 
without change from the GUSTMDL program. In GUSTMDL the coefficients in the differential and 
difference equations were computed one time in the INITIAL BLOCK, since aircraft airspeed VTOT was 
kept constant. In the HARV simulation some of these coefficients must be computed in the DERIVATIVE 
or DISCRETE BLOCKS to accommodate changes in airplane airspeed during the simulated flight. 

In this section the entire HARV simulation code will not be presented for to do so would significantly 
increase the length of the report while adding very little information relative to the turbulence models and 
their implementation. Instead, portions of code which show all of the modifications made to the simulation 
to install the turbulence models will be included. In the code presentation below, to avoid confusion sections 
of simulation code are separated from portions of report text by a line of #’s. 

In the MACRO SECTION MACRO DRMS was inserted after the existing MACRO ALFDYN. The 
code for MACRO DRMS follows and was used to compute statistics needed for evaluation of the performance 
of the three turbulence models. 

########################### ############# 
i 

I**********************************************************************! 
MACRO DRMS ( RMSX , MNX , X, N ) 

i **********************************************************************, 
i 

! MACRO TO COMPUTE SAMPLE STATISTICS OF RANDOM VARIABLE 
! WTB , JCY 2-24-94 

i 

MACRO REDEFINE SUM, SUMSQ, EPS 

CONSTANT SUM = 0 . , SUMSQ = 0 . 

CONSTANT EPS = l.E-22 
SUM = SUM + X 
SUMSQ = SUMSQ + X*X 
MNX = SUM/N 

RMSX = SQRT (MAX ( (SUMSQ - MNX*SUM) / (MAX (N-l . , EPS) ) , 0 . ) ) 

MACRO EXIT 
MACRO END 

########################### ############# 
In the INITIAL BLOCK constants were initialized. The code showing initial values needed by the 
implementation of the three turbulence models follows. Comments indicating equation numbers refer to the 
previously defined equations in this report. 

########################### ########## ### 

RA2DG = 180. /ACOS (-1. ) 

DG2RA = 1./RA2DG 
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CONSTANT UWG = 0 . 
CONSTANT VWG = 0. 
CONSTANT WWG = 0. 

i 

CONSTANT UWND = 0 . 
CONSTANT VWND = 0 . 
CONSTANT WWND = 0 . 


CONSTANT VWNDIC =0., & 

HDWIC = 0 . , & 

VWNDRT = 0 . , & 

HDWRT = 0 . 

######################################## 
Still in the INITIAL BLOCK definitions of turbulence constants and initial values of turbulence 
variables were placed after the END OF/EULER ANGLE/QUATERNION/DIRECTION COSINE/ TRIM 
UPDATE comment in the aircraft simulation. 

######################################## 

T end OF/EULER ANGLE/QUATERNION/DIRECTION COSINE/ TRIM UPDATE.! 

I**********************************************************************! 

i 

,**********************+*★*************************+******+******+*****1 
i START IMPLEMENTATION OF TURBULENCE MODEL WTB,JCY 9-10-97! 

i**********************************************************************! 


CONSTANT 

SDNOIS 

= 1. , 

Sc 


TSAMP 

= .0125, 

Sc 


TURBL 

= 1750 . , 

Sc 


TURBSIG 

= 5. 


CONSTANT 

FILNU 

= 0. , 

Sc 


FILNV 

= 0., 

Sc 


FILNW 

= 0. , 

Sc 


FILNP 

= 0., 

Sc 


BWING 

= 37.4 



CONSTANT FILU 

= 0. , 

Sc 

FILV 

= 0., 

Sc 

FILW 

= 0. , 

Sc 

FILP 

= 0., 

Sc 

FILQ 

= 0. , 

Sc 

FILR 

= 0. 



CONSTANT MILUK = 0 . , & 

MILVK = 0 . , & 

MILWK = 0 . , & 

MILPK = 0 . , & 

MILQK = 0 . , & 

MILRK = 0 . 

I 

######################################## 
Character variable TURBFLG was defined to allow a choice of which turbulence model to be 
implemented during execution of the program. Values of this variable can be ‘CON’, ‘TUS\ or ‘MIL’. 

######################################## 
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CHARACTER TURBFLG*3 
CONSTANT TURBFLG = ' TUS ’ 

CONSTANT EPSLON = l.e-22 

CONSTANT SAMPS = 0. 

i 

LOGICAL PQRFLG 

t 

CONSTANT PQRFLG = .TRUE. 

########################### ############# 
Logical variable PQRFLG was defined to allow a choice to use p-, q-, or r- gusts to turbulence model 
outputs. Values of .True, or .False, should be used for this variable. The following variables were defined 
before the end statement of the INITIAL BLOCK as in the GUSTMDL program. These variables are 
constant and do not depend on velocity. 

########################### ### ########## 

j 

pi = acos ( - 1 . ) 

i 

! irseed should be pos, odd, integer , possibly 8 digits 

i 

! multiply defined irseed CONSTANT IRSEED = 28545269 

GAUSI ( IRSEED ) 

i 

CONSTANT TURBUO = 0 . , & 

TURBXVDO = 0 . , & 

TURBXVO = 0 . , & 

TURBXWDO = 0 . , & 

TURBXWO = 0 . 

i 

CONSTANT TURBRO = 0 . , & 

TURBQO = 0 . , & 

TURBPO = 0 . 

; 

SIGU = TURBSIG ! same as sigv,sigw 

SIGV = TURBSIG 
SIGW = TURBSIG 

SIGP = 1.9 / SQRT (TURBL * BWING) * SIGW ! eq.8 

i 

! APPROXIMATIONS OF EXPECTED VALUES OF STD. DEV. (P AND R COMPONENTS) 

i 

SIGQ = SQRT (pi/ (2 . *TURBL * BWING)) * SIGW ! eq.14 

SIGR = SQRT (2 . *pi/ (3 . *TURBL * BWING)) * SIGW ! eq.17 

; 

J *********★************★*★★****************************★***************! 

END $ ! OF INITIAL ! 

J**********************************************************************I 
I ********************************************************************** | 

I ********************************************************************** \ 

DYNAMIC 

i **★★★★***★*★★ + ★★*******★*****★*******★★★**★*★★★*★★*★■*•**★★★★★* + *■*■****** i 

DERIVATIVE 
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######################################## 
Variables implementing the continuous models were placed in the DERIVATIVE BLOCK of the ACSL 
simulation. They were inserted with the steady-state wind model after the END statement for the linearization 
procedural and before the table look-ups for the HARV aerodynamic model. This placement of the model 
allows the latest turbulence velocities to be included in the calculation of aerodynamic forces and moments. 

######################################## 

i 

NOLLN . . CONTINUE 

; 

END$ ! ENGINE LINEARIZATION PROCEDURAL! 


i COMPUTE BODY FRAME COMPONENTS OF WIND! 

VELWND = VWNDIC + VWNDRT* ( H - 15000 .) /1000 . 

HEADWND = HDWIC + HDWRT* ( H - 15000 .) /1000 . 

XDWND = VELWND * COS ( HEADWND* DG2RA) 

YDWND = -VELWND * SIN ( HEADWND *DG2RA) 

ZDWND = 0 . 

XDDWND = VWNDRT* (HD/1000 . ) *COS ( HEADWND *DG2RA) & 

- VELWND *S IN ( HEADWND *DG2RA) * HDWRT* (HD/1000. ) *DG2RA 
YDDWND = - VWNDRT* (HD/1000 .) *SIN ( HEADWND *DG2RA) & 

-VELWND*COS ( HEADWND *DG2RA) * HDWRT* (HD/1000 . ) *DG2RA 
ZDDWND = 0. 

UWND = CXX* XDWND + CXY* YDWND + CXZ*ZDWND 

VWND = CYX* XDWND + CYY* YDWND + CYZ*ZDWND 

WWND = CZX*XDWND + CZY* YDWND + CZZ*ZDWND 


**★*★★*★** 


Turbulence or Gust components 


UAM 

= U - UWND 

VAM 

= V - VWND 

WAM 

= W - WWND 

UA2M 

= UAM* UAM 

VA2M 

= VAM*VAM 

WA2M 

= WAM* WAM 


######################################## 
The variable VTMINUS defines a “pseudo steady-state” value of true airspeed for use in calculating 
turbulence model coefficients. VTMINUS is the total air speed minus, or without, turbulence. 

######################################## 

VTMINUS = SQRT ( UA2M + VA2M + WA2M ) 

VTOT = VTMINUS 

I 

! use turbl= 1750. for Lu, Lv, Lw 
! use tau for tauu, tauv, tauvw 

i 

TAU = TURBL/VTOT 
TAUQ = 4.0 * BWING /(PI * VTOT) 

TAUR = 3.0 * BWING /(PI * VTOT) 

i 


eq. 3 , 4 , 6 
eq. 13 
eq . 16 
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LP = SQRT (TURBL * BWING) / 2.6 
TAUP = LP/VTOT 


! eg. 10 
! eg. 13 


TURBOMEG A = VTOT/ TURBL 
WP = VTOT/LP 


same as 1 . /TAU, 1 . /TAUV, 1 . /TAUW 
same as 1 . /TAUP 


CONSTANTS FOR CALCULATIONS OF CONTINUOUS W AND V 


! eg. 1,2 


TURBK2U 

K2P 

TURBK1VW 

TURBK2VW 

TURBK3VW 

TURBK4VW 


TURBSIG*SQRT (2 . *TURBOMEGA/TSAMP) 
SIGP * SQRT (2. * WP/TSAMP) 

2 . *TURBOMEGA 
TURBOMEGA*TURBOMEGA 
TAU* SQRT ( 3 . ) 

TURBS I G * SQRT ( TURBOMEGA* * 3 / TSAMP ) 






TURBULENCE MODEL WTB, JCY 2-24-94 
*************** ************************* ****************************** , 
U - COMPONENT 

*****************★*************************************★******+*******, 
used in eq.5 


TURBUD = - TURBOMEGA * TURBU + TURBK2U * FILNU 
TURBU = INTVC (TURBUD , TURBUO) 

********************************************************************** f 
V - COMPONENT 






used in eq.2 


TURBXVDD = - TURBK1VW * TURBXVD - TURBK2VW * TURBXV & 

+ TURBK4VW * FILNV 

TURBXVD = INTVC (TURBXVDD, TURBXVDO) 

TURBXVD I = TURBXVD 

TURBXV = INTVC (TURBXVDI, TURBXVO ) 

TURBV = TURBXV + TURBK3VW * TURBXVDI 

TURBVD = TURBXVDI + TURBK3VW * TURBXVDD 

**********************************************************************! 
W - COMPONENT 




* * * * * i 


used in eg. 2 


TURBXWDD = - TURBK1VW * TURBXWD - TURBK2VW * TURBXW & 

+ TURBK4VW * FILNW 

TURBXWD = INTVC (TURBXWDD, TURBXWDO) 

TURBXWD I = TURBXWD 

TURBXW = INTVC (TURBXWD I, TURBXWO ) 

TURBW = TURBXW + TURBK3VW * TURBXWDI 

TURBWD = TURBXWDI + TURBK3VW * TURBXWDD 

i **********************************************************************! 

! P, Q, R COMPONENTS ADDED 5-15-95 JCY 
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I**********************************************************************! 

! P - COMPONENT 

I *★*****★★*★★★★★★**+*+***★**+***+******★*★****★*****★★★**+★*★*★*★★★*★**! 

I 

! used in eq.7 

i 

TURBPD = - WP * TURBP + K2P * FILNP 
TURBP = I NT VC (TURBPD, TURBPO) 

i ********************************************************************** \^ 

l Q - COMPONENT 

T ********************************************************************** I 
] 

! used in eq.12 

; 

TURBQD = - TURBQ / TAUQ + TURBWD / (TAUQ * VTOT) 

TURBQ = I NT VC (TURBQD, TURBQ 0 ) 

| ********************************************************************** I 

! R - COMPONENT 

I************************************************-*-*********************! 

1 

! used in eq . 15 

t 

TURBRD = - TURBR / TAUR + TURBVD / (TAUR * VTOT) 

TURBR = INTVC (TURBRD, TURBR 0 ) 

I *******★*********★*★**************+****************************★******! 

########################### ############# 
The following code selects the continuous, Tustin, or MIL STD models for use in the aircraft 
simulation. 

########################### ############# 

PROCEDURAL (wwg=wwnd, TURBU, TURBV , TURBW # TURBP, TURBQ, TURBR) 

IF (TURBFLG . EQ . ' TUS 1 .AND. T . GT . 0.) THEN 

UWG = FILU 
VWG = FILV 
WWG = FILW 
PWG = FILP 
QWG = FILQ 
RWG = FILR 

IF (.not. PQRFLG) THEN 
PWG = 0. 

QWG = 0. 

RWG = 0 . 

END IF 

UDWG = FILUD2 
VDWG = FILVD2 
WDWG = FILWD2 

ELSE IF (TURBFLG . EQ . 'MIL' .AND. T .GT. 0.) THEN 
UWG = MILUK 
VWG = MILVK 
WWG = MILWK 
PWG = MILPK 
QWG = MILQK 
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RWG = MILRK 

IF (.not. PQRFLG) THEN 
PWG = 0. 

QWG = 0. 

RWG = 0. 

END IF 

UDWG = FILUD2 
VDWG = FILVD2 
WDWG = FILWD2 

ELSE IF (TURBFLG . EQ . 'CON' .AND. T .GT. 0.) THEN 
UWG = TURBU 
VWG = TURBV 
WWG = TURBW 
PWG = TURBP 
QWG = TURBQ 
RWG = TURBR 

IF (.not. PQRFLG) THEN 


PWG 

= 0 . 

QWG 

= 0 . 

RWG 

= 0 . 

END IF 


UDWG = 

TURBUD 

VDWG = 

TURBVD 

WDWG = 

TURBWD 


ELSE 


UWG = 

0 . 

VWG = 

0 . 

WWG = 

0 . 

PWG = 

0 . 

QWG - 

0 . 

RWG = 

0 . 

UDWG = 

0 

VDWG = 

0 

WDWG = 

0 


END IF 

END $ ] OF PROCEDURAL TO CALCULATE turbulence! 


Compute time rate of change of body frame components of wind. 
Vb(ody) = [L] Ve(arth) 

d/dt Vb = (d/dt [L] ) Ve + [L] (d/dt Ve) 

d/dt [L] = - [omega] [L] Etkin, eq (5.2,11) 
d/dt Vb = - [omega] [L]Ve + [L] (d/dt Ve) 
d/dt Vb = - [omega] Vb + [L] (d/dt Ve) 


UDWND = r *vwg - q*wwg + CXX*XDDWND + CXY*YDDWND + CXZ*ZDDWND 

VDWND = -r*uwg + p*wwg + CYX*XDDWND + CYY*YDDWND + CYZ*ZDDWND 

WDWND = q*uwg - p*vwg + CZX*XDDWND + CZY*YDDWND + CZZ*ZDDWND 

UA = U - UWND - UWG 

VA = V - VWND - VWG 

WA = W - WWND - WWG 

UDA = UD - UDWND - UDWG 

VDA = VD - VDWND - VDWG 

WDA = WD - WDWND - WDWG 
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UA2 = UA*UA 
VA2 = VA*VA 
WA2 = WA*WA 

######################################## 
The variables PA, QA, and RA define instantaneous attitude rates relative to the turbulent atmosphere 

for use in the aerodynamic model table look-ups. 

######################################## 


Calc PA, QA, RA 

PA = P + PWG 
QA = Q + QWG 
RA = R + RWG 

PROCEDURAL (TMPK=H) 

CALL ATMAT62 (TMPK, H, 3 ) 

TMPKOLD = TMPK 
TMPK = TMPK + DELTK 

END $ ! OF PROCEDURAL TO CALCULATE TEMPERATURE IN KELVIN! 
COMPUTE WIND AXES VARIABLES ! 


CSALF = COS(ALF) 
CSBET = COS (BET) 
SNALF = SIN(ALF) 
SNBET = SIN (BET) 


PW = ( P* CSALF* CSBET 

QW = ( -P*CSALF*SNBET 

RW = ( -P* SNALF 

PWDG = PW * RA2DG 
QWDG = QW * RA2DG 
RWDG = RW * RA2DG 
AYW = VT * RW 
AZW = -VT * QW 

######################################## 
The variables PA, QA, and RA were used as parameters in the calls to subroutines SFAERRF and 
AEROINC , which define the H ARV aerodynamic model, to replace the inertial variables P, Q, and R . These 
variables (PA, QA, and RA) were also used in the PROCEDURAL surrounding the calls to subroutines 
SFAERRF and AEROINC. Existing variables ALDGRF, BEDGRF, MACHRF were used in the calls to 
the aero models, but these variables now include the effects of turbulence. 

######################################## 


+ R* SNALF*CSBET + 
- R* SNALF* SNBET + 
+ R* CSALF + 


(Q-ALFD) * SNBET ) 
(Q-ALFD) *CSBET ) 
BETD ) 


PROCEDURAL ( 

FVRA 

, CXRFSF , CYRFSF 
, CDRFSF , CLRFSF 
, cno , cnbeta 

, cnp , clo 


& 

& 

, CZRFSF , C1RFSF , CMRFSF , CNRFSF & 

, TV J P AVG , TV JY AVG , TV JRAVG & 

, cnlef , cntef , cnail , cnrud , cnraero& 
, clbeta , cllef ,cltef ,clail ,clrud & 
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, clr 

/ clp 

, dclf lx 





Sc 

=ALDGRF 

, BEDGRF 

, MACHRF 

,PA 

/ QA 

, RA 

, HRF 

Sc 

, HGC 

, QBAR 

,CWRF 

, BWRF 

, VTRF 



Sc 

, DSR 

, DSL 

, DAR 

, DAL 

, DRR 

, DRL 

,DFR 

Sc 

, DFL 

,DNR 

, DNL 

, DSB 

,DLG 



5c 

, LDBA 

, LQSE 

, LRTE 

, SWRF 

t 



5c 

FGAERO 

,H) 

, TVJANP 

, TV J ANY 

, XSTRAKE 

, DG2RA 

, DNSL 

, DNSR 

5c 


CALL SFAERRF TO USE THE R/T AERO DATA E FOR RIGID MOTION. 
SFAERRF INTERPRETS LEFT RUDDER AS POS T . E . LEFT . THIS IS 
IN ACCORDANCE WITH MCAIR CONVENTION. 

Altitude input to SFAERRF has been changed from pressure altitude 
to true altitude. This was done to match the implementations in 
the DMS, ACES, and Dryden simulations. 

CALL SFAERRF ( & 


CDRFSF 

, CYRFSF 

, CLRFSF 

, C1RFSF 

, CMRFSF 

, CNRFSF 


& 

, ALDGRF 

, BEDGRF 

, MACHRF 

, PA 

/ QA 

, RA 

t H 

5c 

, HGC 

, QBAR 

, CWRF 

, BWRF 

, VTRF 



5c 

, DSR 

, DSL 

, DAR 

, DAL 

r DRR 

, DRL 

, DFR 

5c 

, DFL 

, DNR 

, DNL 

, DSB 

, DLG 



5c 

, LDBA 

, LQSE 

, LRTE 

, cno 

, cnbeta 

, cnlef 

, cntef 

5c 

, cnail 

, cnrud 

, cnraero 

r cnp 

, clo 

, clbeta 

, cllef 

5c 

, cltef 

, clail 

, clrud 

, clr 

, clp 

, dclf lx) 




i 

CALL AEROINC ( 


CDRFSF 

, CYRFSF 

, CLRFSF 

, C1RFSF , 

CMRFSF , CNRFSF, & 

ALDGRF 

, BEDGRF 

, MACHRF 

, 

& 

QBAR 

; SWRF 

, FGAERO 

, TVJANP , 

TV J ANY , & 

DSR 

, DSL 

, DRR 

,DRL 

TVJPAVG, TVJYAVG, & 

TVJRAVG , PA 

/ QA 

, RA 

CWRF , BWRF , & 

VTRF 

, XSTRAKE, DG2RA 

, DNSL 

DNSR , & 

DCNFS 

) 




# # # # 

# # # # # 

# # # # 

# # # # # 

##### ############# 


The Tustin model and the MIL STD model were implemented in the DISCRETE BLOCK GUST which 
computes the turbulence outputs from these models at intervals of 0.0125 seconds (80 Hz). 

########################### ############# 



DISCRETE GUST 



INTERVAL TSGAUS = 0.0125 



GAUSSIAN RANDOM PROCESS 
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IF ( SDNOIS. LE.EPSLON) GOTO LETAl 

FILNU = GAUSS i 

:o. 

, SDNOIS) 

FILNV = GAUSS 

(0. 

, SDNOIS) 

FILNW = GAUSS 

(0. 

, SDNOIS) 

FILNP = GAUSS 

(0. 

, SDNOIS) 


GOTO LETA2 
LETAl . . CONTINUE 
FILNU = 0 . 
FILNV = 0. 
FILNW = 0. 
FILNP = 0. 
LETA2 . .CONTINUE 


************************* 


r************ 




CONSTANTS FOR CALCULATIONS of TURBULENCE VIA TUSTIN TRANSFORM 
TWO P I OVTNU = SQRT(2 . *PI/TSAMP) 


********************************************************************** ! 

U-COMPONENT 

***************************************************** ***************** ! 

constants used in eq.18 

KFIL = TURBSIG* SQRT ( 1 . / (PI * TURBOMEGA) ) 

CFIL = TURBOMEGA/ TAN (TURBOMEGA*TSAMP/2 . ) ! eq.19 

FILK1 = (TURBOMEGA - CFIL) / (TURBOMEGA + CFIL) 

FILK2 = KFIL /(l. + (CFIL/TURBOMEGA) ) 

FILK3 = SQRT(PI * (TURBOMEGA + CFIL ) ) /SDNOIS 
,********************************************************************** ! 
V-COMPONENT and W-COMPONENT 

********************************************************************** ! 

constants used in eq.20 

KVW = SQRT ( 3 . * TURBOMEGA/ ( 2 . * PI ) ) 

WNPC = TURBOMEGA + CFIL 

FILK4 = 2 .* (TURBOMEG A* TURBOMEG A - CFIL*CFIL) / (WNPC*WNPC) 

FILK5 = (TURBOMEGA - CFIL) * (TURBOMEGA - CFIL) / (WNPC*WNPC) 

FILK6 = CFIL + TURBOMEGA /SQRT ( 3 . ) 

FILK7 = 2 . *TURBOMEGA/SQRT ( 3 . ) 

FILK8 = TURBOMEGA/ SQRT ( 3 . ) - CFIL 

FILK9 = KVW* TURBSIG / (WNPC*WNPC) 


used in eq.90 


KFILD1 

potential 


= TURBSIG*SQRT (2 . *TAU/TSAMP) 
alternate (filwd2, filvd2) 
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KWD1 = (1. - TAU*CFIL) / ( 1 . + TAU*CFIL) 

KWD2 = TURBSIG*SQRT (TAU/TSAMP) *CFIL/ ( 1 . +TAU*CFIL) **2 
CTWSR3 = SQRT ( 3 . ) *TAU*CFIL 


P- COMPONENT 


used in eq.21 

PFIL = SIGP* SQRT ( 2 . / (TSAMP * WP) ) 
PCFIL = WP / TAN ( WP *TSAMP/2 . ) 

PK1 = (WP - PCFIL) / (WP + PCFIL) 
PK2 = PFIL / ( 1 . + ( PCFIL/WP ) ) 


! eq.22 


Q-COMPONENT and R-COMPONENT 


used in eq.23 


CFILQ = (1. / TAUQ) / TAN (TSAMP / (2. * TAUQ) ) ! eq . 24 

QK1 = (1. - CFILQ * TAUQ) / (1. + CFILQ * TAUQ) 

QK2 = (CFILQ / VTOT) / (1. + CFILQ * TAUQ) 

used in eq.25 

CFILR = (1. / TAUR) / TAN (TSAMP / (2. * TAUR) ) ! eq. 26 

RK1 = (1. - CFILR * TAUR) / (1. + CFILR * TAUR) 

RK2 = (CFILR / VTOT) / (1. + CFILR * TAUR) 






constants for MIL STD calcs 8-28-97 

used in eqs . 30 - 35 


tovtau = tsamp/tau 
tovtau2= 2.* tovtau 
tovtau4= 4.* tovtau 

! (for 

u, w, v) 

tovtaup = tsamp/taup 
tovtaup2= 2.* tovtaup 

! (for 

p) 

tovtauq = tsamp/tauq 

I ( for 

q) 

tovtaur = tsamp/taur 

i ( for 

q) 

milkl = (1. - tovtau) 
milk2 = (1. - tovtau2) 
milk3 = sqrt ( tovtau2 ) 
milk4 = sqrt ( tovtau4 ) 

! ( for 

u, w, v) 

milk5 = (1. - tovtaup) 
milk6 = sqrt ( tovtaup2 ) 

! ( for 

p) 

milk7 = (1. - tovtauq) 

I ( for 

q) 

irtilk8 = (1. - tovtaur) 
piov4b=pi / (4 . *bwing) 
piov3b=pi / (3 . *bwing) 

! ( for 

r) 
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**★***★*★***★************************'* 

TURBULENCE VIA BILINEAR TRANSFORM, OR TUSTIN TRANSFORM 


********************************************* 




U- COMPONENT 

**********************************************************************! 


UFILK = FILNU 
IF (T . EQ. 0.) THEN 
UFILKM1 = 0. 
GUFILKM1 = 0. 
FILUDKMl = 0. 
FILUD2KM1 = 0. 

END IF 


eq . 18 

GUFILK = -FILK1*GUFILKM1 + TWOPIOVTNU*FILK2* (UFILK + UFILKM1 ) 




r ★ ★ ★ ★ ★ 1 


■ fr ****************'****************** 1 

DIGITAL IMPLEMENTATION OF DERIVATIVE CALCULATIONS OF U-COMPONENT 
*********************************************************************** 


eq . 27 

FILUD = (GUFILK - GUFILKMl ) /TSAMP 

eq . 2 8 

FILUD2 = - KWDl * f ilud2kml 

+ (kfildl*CFIL/ (1. + tau*CFIL) ) 

*(ufilk - ufilkml) 


FILUDKMl = FILUD 
FILUD2KM1 = FILUD2 
UFILKM1 = UFILK 
GUFILKMl = GUFILK 
FILU = GUFILK 


************************** 


*************************** 




********* 


I 


V-COMPONENT ****************1 


VFILK = FILNV 
IF (T .EQ. 0.) THEN 
VFILKM1 = 0. 
VFILKM2 = 0 . 
GVFILKMl = 0. 
GVFILKM2 = 0 . 
FILVDKM1 = 0. 
FILVD2KM1 = 0. 
FILVD2KM2 = 0. 

END IF 

FILVKMl = GVFILKMl 


eq . 20 
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GVFILK = - FILK4*GVFILKM1 - FILK5 *GVFILKM2 & 

+ TW0PI0VTNU*FILK9* (FILK6*VFILK + FILK7 *VFILKM1 & 

+ FILK8*VFILKM2) 

**********************************i,i,*i,i'* i ' i ,i i ,* i 'i ri , irici , iri ' i ' i ' i ' i ' iril . i!irir i ri(i '** i ' i ' i ' | 

DIGITAL IMPLEMENTATION OF DERIVATIVE CALCULATIONS OF V-COMPONENT 
*****************4****************H**tHt*********t***i*H***********| 


eq . 27 


FILVD = (GVFILK - GVFILKM1) /TSAMP 


eq. 2 9 


FILVD2 = -2 . *KWD1*FILVD2KM1 & 

- KWD1 * *2 *FILVD2KM2 & 

+ kwd2 * ( ( 1 . +CTWSR3 ) *VFILK & 

- (2 . *CTWSR3 ) *VFILKM1 - ( 1 . -CTWSR3 ) *VFILKM2 ) ! 5-12-97 

i 

GVFILKM2 = GVFILKM1 
GVFILKM1 = GVFILK 
VFILKM2 = VFILKM1 
VFILKM1 = VFILK 
FILVD2KM2 = FILVD2KM1 
FILVD2KM1 = FILVD 2 
FILVDKM1 = FILVD 
FILV = GVFILK 

! W-COMPONENT 

i***********t**n**mtt*i*******H***4******tt**** t ** 4 * tt * JHtHnlHni | 

WFILK = FILNW 
IF (T .EQ. 0.) THEN 
WFILKM1 = 0. 

WFILKM2 = 0. 

GWFILKM1 = 0 . 

GWFILKM2 = 0. 

FILWDKM1 = 0. 

FILWD2KM1 = 0. 

FILWD2KM2 = 0. 

END IF 

FILWKM1 = GWFILKM1 

! eq.20 

GWFILK = - FILK4 *GWFILKM1 - FILK5 *GWFILKM2 & 

+ TWOPIOVTNU*FILK9* (FILK6*WFILK + FILK7*WFILKM1 & 

+ FILK8*WFILKM2 ) 

DIGITAL IMPLEMENTATION OF DERIVATIVE CALCULATIONS OF W-COMPONENT 
* ***********************************************************^*^^^^^^^^1 

eq. 27 


FILWD = (GWFILK - GWFILKM1 ) /TSAMP 
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eq. 29 


FILWD2 = -2 . *KWD1*FILWD2KM1 & 

- KWD1**2*FILWD2KM2 & 

+ kwd2*( (1.+CTWSR3) *WFILK & 

- (2 . *CTWSR3) *WFILKM1 - ( 1 . -CTWSR3 ) *WFILKM2 ) ! 5-12-97 

GWFILKM2 = GWFILKM1 
GWFILKMl = GWFILK 
WFILKM2 = WFILKM1 
WFILKM1 = WFILK 
FILWD2KM2 = FILWD2KM1 
FILWD2KM1 = FILWD2 
FILWDKMl = FILWD 
FILW = GWFILK 

*****★***** + ******★******** + ******************************•*** + ******** t 

P, Q, R COMPONENT ADDED 5-15-95 JCY 

★ *************** + ******************************************** + ******** | 
P- COMPONENT 

**************************************************+**********+**+*****t 

PFILK = FILNP 
IF (T .EQ. 0.) THEN 
PFILKM1 = 0. 

GPFILKMl = 0. 

END IF 

GPFILK = - PK1 * GPFILKMl + PK2* (PFILK + PFILKM1 ) ! eq . 21 

PFILKM1 = PFILK ! 5-2-97 

GPFILKMl = GPFILK ! 5-2-97 

FILP = GPFILK 

**********************★************* + ***★*************★**************★1 

Q- COMPONENT 

★ ★★a-****************************************************************** I 
IF (T .EQ. 0.) THEN 
FILQKMl = 0. 

END IF 

FILQ = -QK1* FILQKMl + QK2*(FILW - FILWKM1 ) ! eq. 23 

FILQKMl = FILQ 

a********************************************************************* I 

R-COMPONENT 

**********************************************************************i 

IF <T .EQ. 0.) THEN 
FILRKMl = 0. 

ENDIF 

FILR = -RK1* FILRKMl + RK2* (FILV - FILVKM1 ) ! eq . 25 

FILRKMl = FILR 

********************************************************************** i 

ALPHA , BETA COMPONENTS ADDED 5-7-97 

a*********************************************************************! 

FILA = FILW/VTOT ! ALPHA COMPONENT 
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FILB = FILV/VTOT 1 BETA COMPONENT 

FILAD = FILWD/VTOT 
FILBD = FILVD/VTOT 

**★***** + ★*** + **★★★*★★★*********★**★**★*******★★★★*★****★★***★★**★■*■*** t 

ALPHA , BETA COMPONENTS ALTERNATES ADDED 6-19*97 

********************************************************************** j 

FI LAD 2 = FILWD2/VTOT 
FILBD2 = FILVD2/VTOT 


*******************************^**************************************| 
MIL STD IMPLEMENTATION OF CALCS OF u,v,w,p,q,r- COMPONENT 

8-28-97 

********** + ***★**★***★************★* + *********************★******★*★** j 

IF (T .EQ. 0.) THEN 
MILUKM1 = 0. 

MILVKM1 = 0. 

MILWKM1 = 0. 

MILPKM1 = 0. 

MILQKM1 = 0. 

MILRKM1 = 0. 

END IF 


MILUK = milkl*MILUKMl + sigu*milk3*ufilk i eq.30 
MILVK = milk2*MILVKMl + sigv*milk4 *vf ilk ! eq.31 
MILWK = milk2*MILWKMl + sigw*milk4*wf ilk ! eq.32 
MILPK = milk5*MILPKMl + sigp*milk6*pf ilk ! eq.33 
MILQK = milk7*MILQKMl + piov4b* (milwk-milwkml ) ! eq.34 
MILRK = milk8 *MILRKM1 + piov3b* (milvk-milvkml ) ! eq.35 

Calc derivatives 9-23-97 


MILUDK = (MILUK - MILUKM1 ) /TSAMP 
MILVDK = (MILVK - MILVKM1 ) /TSAMP 
MILWDK = (MILWK - MILWKM1 ) /TSAMP 

! save past values 

MILUKM1 = MILUK 
MILVKM1 = MILVK 
MILWKM1 = MILWK 
MILPKM1 = MILPK 
MILQKMl = MILQK 
MILRKM1 = MILRK 

i 

J ********************************************************************** 1 

END ! GUST DISCRETE 

i *★**★*★**★★**★*★★★****★***★**★****★***★★*★★*******★******★**** + ******★ i 
i **★★*★★*★★★★*★**★*********★★*★*★★**★★★★★★*★★★* + **★★■*■*■*•***** + ★★*★★*★★** i 

DISCRETE DISCRMS 
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■a.***.*.**.*..*.******************** 




DISCRETE TO COMPUTE RMS OF TURBULENCE 

INTERVAL TRMS = 0.0125 
SAMPS = SAMPS + 1. 

DRMS (DRMSUWG, DMNUWG, UWG , SAMPS) 
DRMS ( DRMSVWG , DMNVWG , VWG , SAMPS ) 
DRMS (DRMSWWG, DMNWWG, WWG , SAMPS) 

DRMS (DRMSPWG, DMNPWG, PWG , SAMPS) 
DRMS ( DRMSQWG , DMNQWG , QWG , SAMPS) 
DRMS ( DRMSRWG , DMNRWG, RWG , SAMPS) 

DRMS (DRMSUDWG, DMNUDWG, UDWG , SAMPS) 
DRMS ( DRMSVDWG , DMNVDWG , VDWG , SAMPS ) 
DRMS (DRMSWDWG, DMNWDWG, WDWG , SAMPS) 


MIL STD calcs, u,v,w,p,q,r extra if needed 

DRMS (DRMSMILUK, DMNMILUK, MILUK, SAMPS) 
DRMS (DRMSMILVK, DMNMILVK, MILVK, SAMPS) 
DRMS (DRMSMILWK, DMNMILWK, MILWK, SAMPS) 
DRMS (DRMSMILPK, DMNMILPK, MILPK, SAMPS) 
DRMS (DRMSMILQK, DMNMILQK, MILQK, SAMPS) 
DRMS (DRMSMILRK, DMNMILRK, MILRK, SAMPS) 


END ! of DISCRMS DISCRETE 

*********************************************************************! 


########################### ############# 


Time History Plots 

Numerous runs were made with the modified HARV aircraft simulation to evaluate the performance of 
the implemented turbulence models. The 20-second runs were made with five trim cases. Some runs set the 
PQRFLG variable to .True, to use effects of p, q, and r gust effects, while other runs had a value of .False, 
for the PQRFLG variable. All runs used the same random seed for noise generation. Time history plots of 
comparisons of the continuous, Tustin, and MIL STD models for each trim case are presented below in 
figures 15 through 59. Attitude and air data variables are also plotted. The five trim cases plotted were for a 
values of 5°, 25°, 35°, 55°, and 60°. Figures 15 through 23 are the time histories for the a = 5° trim case. 
Figures 24 through 32 are the time histories for the a = 25° trim case. Figures 33 through 41 are the time 
histories for the a = 35° trim case. Figures 42 through 50 are the time histories for the a = 55° trim case. 
Figures 51 through 59 are the time histories for the cc = 60° trim case. 
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Figure 19. Comparison of continuous and Tustin air data and accelerations for a = 5°, seed no. 1 
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(c).- Yaw rate. 
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Figure 20. Comparison of Tustin attitude rates and a 


phidg, deg <1 d g- deg/sec 







ayft, ft/sec 2 vt » ft^ sec alfdg, deg 




















ayft, ft/sec 2 vt, ft/sec alfdg, deg 



5 10 

Time, sec 

(a).- Angle of attack. 


5 10 15 

Time, sec 

(b).- Sideslip angle. 



10 

Time, sec 

(c).- True airspeed. 



5 10 1 

Time, sec 

(d).- X-axis acceleration. 



5 10 15 

Time, sec 

(e).- Y-axis acceleration. 



MIL STD - w/ PQR 

MIL STD - w/o PQR 




10 

Time, sec 


(e).- Y-axis acceleration. (f).- Z-axis acceleration. 

Figure 23. Comparison of MIL STD air data and accelerations w/ and w/o PQR gusts for a = 5°, seed no. 1 
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Figure 24. Comparison of continuous, Tustin, and MIL STD model turbulence for a = 25°, seed no. 1. 
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Figure 26. Comparison of continuous and Tustin air data and accelerations for a = 25°, seed no. 1 . 


1 


I 













vt, ft/sec 



109 














1 


-Tustin - w / PQR I 
Tustin - w/o PQRl 


o a. 




10 

Time, sec 


Tustin - w/ PQR 
Tustin - w/o PQR 



(a),- Roll rate. 


10 

Time, sec 
(b).- Pitch rate. 



M - Tustin - w/ PQR ( 
Tustin - w/o PQRl 


10 

Time, sec 
(c).- Yaw rate. 



— . — . — . — - — i — - — • — 1 — « — ! — - — ' — ' — 


Tustin - w/ PQR 


Tustin - w/o PQR 


j \ 


j K i 


J | \ ; 


1 " i 

! i . . i 

i — * — i — i — i i . i ' 


5 10 15 

Time, sec 

(e).- Pitch angle. 


0.4 

0.2 

0 

OD 

•o -0.2 
d) 

i -0.4 
- 0.6 


10 

Time, sec 


(d).- Bank angle. 


"O’® “ I Tustin - w/ PQR I 

’ I Tustin - w/o PQRl 


10 

Time, sec 



(e).- Pitch angle. (f).- Heading. 

Figure 29. Comparison of Tustin attitude rates and angles w/ and w/o PQR gusts for a = 25 °, seed no. 1 . 
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Figure 30. Comparison of Tustin air data and accelerations w/ and w/o PQR gusts for a = 25°, seed no. 1. 
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Figure 34. Comparison of continuous and Tustin attitude rates and angles for a = 35 , seed no. 1. 
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Figure 35. Comparison of continuous and Tustin air data and accelerations for oc = 35° seed no. 1. 
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Comparison of continuous and MIL STD attitude rates and angles for a = 35 , seed no. 1 . 
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Figure 37. Comparison of continuous and Tustin air data and accelerations for a = 35° seed no 1 
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Figure 38. Comparison of Tustin attitude rates and angles w/ and w/o PQR gusts for a = 35°, seed no. 1 . 
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Figure 39. Comparison of Tustin air data and accelerations w/ and w/o PQR gusts for a = 35°, seed no. 1 
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Figure 42. Comparison of continuous, Tustin, and MIL STD model turbulence for a = 55°, seed no. 1 
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Figure 43. Comparison of continuous and Tustin attitude rates and angles for a = 55°, seed no. 1 
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Figure 46. Comparison of continuous and Tustin air data and accelerations for a = 55°, seed no. 1 . 
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Figure 47. Comparison of Tustin attitude rates and angles w/ and w/o PQR gusts for a = 55°, seed no. 1. 
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Figure 51. Comparison of continuous, Tustin, and MIL STD model turbulence for a = 60°, seed no. 1 
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Figure 52. Comparison of continuous and Tustin attitude rates and angles for a = 60°, seed no. 1. 
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Figure 53. Comparison of continuous and Tustin air data and accelerations for a = 60°, seed no. 1 . 
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Figure 59. Comparison of MIL STD air data and accelerations w/ and w/o PQR gusts for a = 60°, seed no. 1 
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